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FOREWORD 


This  is  the  final  technical  report  submitted  by  Barron  Associates,  Inc., 
Stanardsville,  Virginia,  under  Phase  I  SBIR  Contract  DAA21-87-C-0140  ("Use 
of  Polynomial  Networks  to  Improve  on  Control  Efficiency  of  Maneuverable 
Projectile")  with  the  U.  S.  Army  Armament,  Munitions  and  Chemical 
Command,  Picatinny  Arsenal,  New  Jersey.  The  work  reported  was  performed 
during  the  period  July  29  -  December  28, 1987.  The  Principal  Investigator  was 
Roger  L.  Barron,  who  is  responsible  for  the  optimum-path-to-go  guidance 
concept.  The  programming  and  numerical  analyses  were  performed 
primarily  by  John  F.  Elder  IV.  Dean  W.  Abbott  assisted  in  the  work. 

The  Contracting  Officer's  Representatives  in  this  project  were  Thomas 
D.  Hutchings  SMCAR-CCS-C  (July  29  -  October  28)  and  Grant  W.  Manning,  Jr., 
SMCAR-CCS-A  (October  29  -  completion).  The  authors  are  indebted  to  them 
and  to  Modesto  J.  Barbarisi,  SMCAR-SCF-RE,  Cliff  Wilkins,  and  A1  Rahe  for 
technical  data,  encouragement,  guidance,  and  support. 
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1. 


SUMMARY  AND  INTRODUCTION 


1,1.  Summary 


An  optimum-path-to-go  (OPTG)  guidance  capability  has  been  designed 
in  preliminary  form  and  investigated  to  demonstrate  its  feasibility  for 
command  adjusted  trajectory  (CAT)  projectiles.  The  principal  design  tools 
used  are  the  calculus  of  variations  and  an  advanced  algorithm  for  synthesis  of 
polynomial  networks  from  simulation  data.  Differential  equations 
governing  optimum  two-point  boundary-value  projectile  guidance  solutions 
for  maneuvering  targets  are  derived  via  the  calculus  and  initialized  in  real¬ 
time  calculations  using  an  adaptively-synthesized  polynomial  network 
(APN).  Application  of  OPTG  guidance  to  120mm  projectiles  has  been 
simulated.  In  particular,  comparisons  have  been  made  of  the  performance  of 
the  OPTG  guidance  when  it  is  blended  with  a  predictive  proportional 
navigation  (baseline)  guidance  law.  In  the  blended  law,  pure  OPTG  guidance 
is  initialized  at  t  =  0  and  used  until  estimated  time-to-go  is  less  than  or  equal 
to  1.0  sec.  (No  re-initialization  of  the  OPTG  guidance  is  used  after  t  =  0.)  The 
baseline  gxiidance  is  then  employed  after  time-to-go  becomes  less  than  1.0  sec. 
until  impact  (or  closest  point  of  approach). 


Block  diagrams  of  the  OPTG  and  baseline  guidance  laws  are  presented  below 
in  Figures  1  and  2.  Figure  l.a  illustrates  an  arrangement  in  which  the  means 
for  intercept  point  prediction  are  separate  from  those  for  computing  the 
initialization  (and  re-initialization)  solutions,  which  is  the  case  investigated 
in  the  present  study.  Figure  l.b  shows  the  means  for  intercept  point 
prediction  combined  with  those  for  the  initialization/re-initialization 
computations,  which  is  the  configuration  recommended  for  future 
development.  Figure  2  (in  Section  5)  presents  a  block  diagram  of  the  baseline 
guidance  law. 


Simulations  of  the  blended  guidance  law  and  pure  OPTG  guidance  have 
demonstrated  the  following  performance  benefits  for  target  trajectories 
having  a  small  level  of  uncertainty  concerning  the  target  accelerations  and 
intercept  coordinates: 
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Figure  l.a:  OPTG  Guidance  Computations 
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Figure  l.b:  Condensed  OPTG  Guidance  Computations 


1.  The  average  closest-point-of-approach  (CPA)  is  2.0  times  less  for 
the  blended  guidance  law  than  for  the  baseline. 

2.  The  area  in  the  miss-distance  plane  reachable  with  a  CPA  less 
than  6  ft.  is  3.5  times  greater  for  the  blended  guidance  law  than 
for  the  baseline. 


3.  For  all  values  of  required  divert,  the  average  improvement 
ratio,  defined  as  CPA  for  the  baseline  divided  by  CPA  for  the 
blended  guidance  law,  is  greater  than  1.0.  For  required  divert 
distances  between  22  and  54  ft.,  the  ratio  is  greater  than  2.0:  and 
between  30  and  51  ft,  the  ratio  is  greater  than  5.0. 

4.  For  the  pure  OTG  guidance,  a  given  divert  is  achieved  with  10 
fewer  squib  firings  on  the  average  than  are  used  with  the 
baseline  guidance.  When  the  blended  guidance  law  is  employed, 
squib  utilization  is  equal  for  the  blended  guidance  law  and  the 
baseline.  (It  is  believed  that  re-initialization  of  the  OPTG 
guidance  after  t  =  0  in  the  blended  guidance  law  will  produce 
most  of  the  squib-usage  reduction  of  pure  OPTG  guidance.) 

1.2  Background 

When  using  conventional  guidance  (such  as  the  baseline  guidance 
law),  a  critical  consideration  is  the  establishment  of  the  threshold  used  for  the 
squib  fire/no-fire  decisions.  If  this  threshold  is  set  too  low,  the  system  may 
use  an  excessive  number  of  squibs  when  the  total  divert  required  is  small  or 
even  negligible.  To  illustrate  this  point,  a  simulation  of  the  baseline  guidance 
has  been  performed  for  a  case  in  which  no  divert  was  required;  that  is,  the 
gun  was  perfectly  aimed  and  purely  ballistic  motion  of  the  projectile  would 
score  an  exact  hit  against  a  non-moving  airborne  target  (such  as  a  helicopter 
in  stationary  hover).  With  the  decision  threshold  of  the  baseline  law  set  to 
zero,  the  projectile  used  almost  all  of  its  40  squibs.  (It  became,  in  effect,  a  self- 
disturbing  as  well  as  self-correcting  system,  a  system  triggered  by  the  slightest 
noise.)  Although  the  miss  distance  was  nearly  zero,  the  zero-level  threshold 


clearly  produced  a  design  that  was  profligate  in  its  use  of  maneuver  resources 
-  a  design  that  could  be  quite  readily  defeated  by  target  maneuvering. 

On  the  other  hand,  if  the  decision  threshold  is  set  too  high,  the  baseline 
guidance  fails  to  react  as  early  in  the  engagement  as  may  be  necessary  if  a  large 
divert  is  required.  As  a  further  illustration,  the  threshold  of  the  baseline 
guidance  was  next  set  to  a  moderate  positive  level  so  as  to  reduce  unnecessary 
squib  firings  for  the  case  just  described  by  about  50  percent.  Miss  distance 
remained  small  as  long  as  little  or  no  divert  was  required,  but  the  design 
exhibited  increased  miss  distances  when  stressed;  that  is,  required  to  approach 
its  maximum  total  divert  capabilities.  The  conclusion:  threshold  selection 
for  conventional  guidance  is  at  best  a  compromise  between  conditions  of 
small  and  large  required  divert.  It  can  be  shown  also  that  this  threshold 
selection  is  a  compromise  between  conditions  of  small  and  large  uncertainties 
(in  measurements  of  system  states  and  predictions  of  target  motions). 

An  important  attribute  of  OPTG  guidance  is  that  the  selection  of  its 
decision  threshold  is  independent  of  the  amount  of  divert  that  will  be 
required.  To  demonstrate  this,  the  OPTG  giudance  was  simulated  for  the  two 
extreme  cases  mentioned  above  (no  divert  required  and  maximum  divert 
required).  The  threshold  was  kept  at  zero  for  the  OPTG  tests.  When  no 
divert  was  required,  no  squibs  were  fired.  When  maximum  divert  was 
required,  the  OPTG  guidance  began  squib  firing  at  the  earliest  available 
opportunity  (as  afforded  by  the  cyclical  variation  of  roll  attitude),  and  kept 
firing  at  each  succeeding  opportunity  until  the  supply  of  squibs  was 
exhausted.  Clearly,  no  strategy  for  guidance  could  have  produced  a  greater 
divert.  Between  the  limits  of  zero  and  maximum  divert,  the  OPTG  guidance 
used  a  number  of  squibs  that  was  approximately  in  proportion  to  the  required 
divert.  At  no  time  was  there  any  question  regarding  the  choice  of  decision 
threshold. 

For  the  OPTG  guidance  system,  the  decision  threshold  can  be  set 
entirely  on  the  basis  of  the  expected  levels  of  noise  and  target  maneuver 
uncertainty.  If  the  threshold  is  placed  above  zero  in  consideration  of  these 
expected  levels,  there  will  be  some  loss  of  maximum  divert  capability.  But 
the  OPTG  guidance  is  appropriately  nonlinear,  and  small  decision  thresholds 


incorporated  to  de-sensitize  it  for  realistic  levels  of  noise  and  uncertainty 
have  an  almost  negligible  influence  on  the  maximum  divert  capability. 
Conversely,  the  baseline  guidance  law  is  proportional  in  its  action,  and  the 
same  threshold  that  adequately  de-sensitizes  both  it  and  the  OPTG  guidance 
to  the  effects  of  noise  and  target  motion  uncertainties  can  be  crippling  to  the 
divert  capability  of  the  baseline  guidance. 

1.3  Optimum-Pa  th-To-Go-Guidance 

Optimum-path-to-go  guidance  provides  optimum,  real-time,  two- 
point  boimdary-value  guidance.  It  periodically  computes  a  new  best  trajectory 
from  the  actual  position  and  velocity  state  of  the  weapon  to  a  designated  final 
position.  In  application  to  CAT  projectiles,  the  OPTG  method  offers  two 
particularly  important  benefits.  First,  by  maneuvering  throughout  the 
engagement  to  establish  the  designated  final  state,  the  magnitude  of  required 
terminal  acceleration  is  minimized,  greatly  reducing  miss  distance.  Second, 
the  OPTG  guidance  realizes  optimum  performance  between  specified 
boundary  conditions  and  thus  enlarges  the  achievable  boundaries  in 
comparison  with  traditional  guidance. 

OPTG  guidance  uses  polynomial  networks,  evolved  from  neural 
networks  beginning  in  the  early  1960s  [1-3,  5-8,  12,  16],  to  implement,  via 
algebraic  inverse  modeling,  time-varying  solutions  for  the  weapon  steering 
commands  as  functions  of  initial,  i.e.,  "time-now",  and  designated  future 
trajectory  states.  (In  trajectory  prediction,  information  about  initial  states  and 
the  steering  commands  is  used  to  compute  a  future  state  or  states.  In  inverse 
modeling,  the  boundary  conditions  are  supplied  and  the  commands  that 
satisfy  those  conditions  are  computed.)  As  do  other  neural  networks,  the 
polynomial  network  learns  from  experience.  However,  the  learning  process 
is  completed  off-line  during  the  design  phase.  The  network  that  is  used  in 
the  weapon  thus  has  fixed  parameters,  and  its  behavior  is  fully  known 
(primarily  from  simulations  but  cilso  from  flight  testing)  before  operational 
usage  occurs. 

The  calculus  of  variations  is  employed  in  design  of  the  guidance  law  so 
as  to  achieve  optimality  and  the  benefits  of  a  compact  representation  of  the 
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optimum  trajectories.  The  procedure  is  first  to  derive  the  governing 
variational  equations,  which  comprise,  chiefly,  a  family  of  first-order, 
nonlinear,  time-varying  adjoint  differential  equations  involving  the 
Lagrange  multipliers  of  the  formulation.  The  derivation  is  based  upon  an 
integral  performance  criterion  having  a  variable  upper  limit  (final  time)  of 
integration.  Next,  an  extensive  data  base  of  optimum  trajectories  is 
computed.  Each  entry  in  this  data  base  consists  of  an  initial  condition  (i.e.,  the 
initial  position  and  velocity  states),  a  set  of  initializing  Lagrange  multiplier 
values  (which  determine  the  specific  optimal  trajectory  flown  as  the  adjoint 
equations  are  integrated),  and  the  resulting  final  condition  (final  states).  The 
data  base  may  be  thought  of  as  a  field  of  extremals  between  large  envelopes  of 
admissible  initial  and  final  condition  pairs.  The  Algorithm  for  Synthesis  of 
Polynomial  Networks  (ASPN)  is  then  used  to  create,  off-line,  a  nonlinear 
polynomial  network  transformation  (called  an  adaptively-synthesized 
polynomial  network  or  APN)  that  maps  observed  initial  and  desired  final 
states  into  initializing  multiplier  values  that  produce  an  optimum  path 
between  the  prescribed  states.  Within  the  weapon,  the  APN  may  be 
interrogated  periodically  to  obtain  optimum-path-to-go  solutions,  and  the 
adjoint  equations  are  integrated  to  compute  steering  commands  until  the 
next  re-initialization  is  performed. 

Potential  further  uses  of  polynomial  networks  are  to  (1)  predict  the 
target  trajectory  and  (2)  estimate  sensitivities  of  weapon  predicted  final  states 
to  changes  in  initial  values  of  the  Lagrange  multipliers,  the  latter  providing 
the  means  for  closed-loop  vernier  adjustments  of  guidance  commands  so  as 
to  null  the  predicted  final  error. 

ASPN  [12]  evolves  the  network  during  the  off-line  design,  beginning 
with  the  simplest  structure  (each  output  connected  directly  to  the  best 
corresponding  input)  and  growing  to  an  appropriately  complex  structure. 
The  final  network  typically  involves  several  layers  of  nodal  algebraic 
elements,  each  having  one,  two,  or  three  inputs  obtained  from  a  previous 
layer  and  providing  a  nonlinear  analytic  function  that  includes  cross-products 
of  its  inputs.  An  information  theoretic  modeling  criterion  is  used  by  ASPN 
during  the  synthesis  process  to  select  the  most  relevant  inputs  and  combine 
them  in  suitable  elements  and  in  a  suitable  overall  network  transformation 


that  is  learned  during  synthesis.  The  criterion  also  regulates  the  allowed 
complexity  of  the  network  to  avoid  overfitting  the  data  (over-specialization). 
The  ASPN  procedure  has  been  undergoing  development  for  27  years  and  has 
been  extensively  proven  in  signal  processing,  prediction,  and  control 
applications. 

1.4  Work  Accomplished 

The  specific  accomplishments  of  this  Phase  I  work  are: 

1.  The  three-degree-of-freedom  (3 DOF)  translational  equations  of 
motion  for  a  120  mm  maneuverable  CAT  projectile  have  been 
written  and  programmed  in  a  simulation  for  trajectory  data  base 
generation.  This  simulation  includes  provisions  for  target 
maneuvers  and  maneuver  uncertainties. 

2.  A  baseline  predictive  guidance  law  has  been  formulated  and  used 
for  comparative  analysis  purposes. 

3.  The  governing  equations  for  optimum-path-to-go  (OPTG)  guidance 
of  the  CAT  projectile  have  been  derived  using  an  extended  calculus 
of  variations  methodology.  Performance  sub-criteria  in  this 
derivation  of  the  two-point  boundary-value  guidance  solution 
include  (a)  minimization  of  thruster  firings,  (b)  maximization  of 
projectile  range,  (c)  maximization  of  projectile  divert  capability,  (d) 
conservation  of  projectile  kinetic  energy,  and  (e)  trajectory  shaping 
factors  associated  with  the  levels  of  uncertainty  regarding  system 
states  and  future  maneuvers  of  the  target. 

4.  The  OPTG  guidance  equations  have  been  incorporated  in  the  3DOF 
simulation  and  numerical  studies  performed,  producing  the 
resuilts  summarized  in  Section  1.1  above. 

5.  Processor  throughput  required  to  implement  the  OPTG  guidance 
has  been  estimated.  (The  throughput  requirement  is 
approximately  5.3  kFLOPS.) 
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CRITERION  OF  PERFORMANCE 


The  Lagrange  formulation  of  the  performance  criterion  involves,  in 
this  application,  the  minimization  of  the  definite  integral 

tf 

I  =  [  Fdt  2:1 

ti 

where  tj  is  the  time  of  guidance  computation  (at  guidance  initiation  t^  =  tQ) 
and  t£  is  the  final  (burst  or  impact)  time.  From  the  mathematical  viewpoint, 
time  tj  is  fixed,  but  t£  is  variable  and  subject  to  optimization.  The  integrand 
(Hamiltonian)  may  be  written  as 

n'  j’ 

F  =  £  W„G„  t  ^  Xjfj  2:2 

n  =  0  j  =  1 

in  which  the  Wj^'s  are  constant  Lagrange  multipliers  that  serve  as  criterion 
weights;  the  G^j’s  are  the  integrands  of  isoperimetric  conditions  that  comprise 
sub-criteria  performance  measures;  and  the  Xj's  are,  in  general,  time-varying 
Lagrange  multipliers  that  introduce  the  anholonomic  (non-integrable) 
equations  of  motion  and  kinematic  relationships. 

It  is  assumed  that  F(t)  is  continuous  and  has  continuous  second  partial 
derivatives  with  respect  to  the  state  variables  of  the  physical  system  to  be 
guided. 

The  following  sub-criteria  are  proposed. 

2.1  Guidance  Effort  Penalty 

Let  N(t)  represent  the  cumulative  number  of  squib  firings  by  the 
projectile  as  of  time  t.  Then,  if  no  squib  firing  occurs  at  t,  N(t)  is  zero,  and  if  a 
squib  is  fired,  N(t)  is  unity  over  the  short  finite  interval  (t,  t  +  At)  between  t 


and  the  next  decision  point.  Note  that  N(t)  is  continuous  even  though  N(t) 
can  exhibit  jump  discontinuities  between  the  levels  zero  and  one. 


The  total  guidance  effort  over  the  flight  of  the  projectile  is  defined  as 
tf 

f  N^dt 
to 


In  the  context  of  optimum-path-to-go  guidance,  the  guidance  effort  to  be 
expended  over  the  time  remaining  as  of  a  solution  time  tj  is 

tf 

f  N2dt 

ti 

This  integral  is  seen  to  be  independent  of  when  squibs  are  fired  and  to  vary 
only  with  the  number  of  firings  over  the  interval  (tj,  tf).  It  may  be 

advantageous  to  fire  the  squibs  early  in  the  engagement,  so  as  to  obtain 
maximum  divert  capability  for  the  expenditure  of  resources.  Conversely,  it 
may  be  advantageous  to  fire  the  squibs  late  in  the  engagement,  so  as  to  use 
resources  when  there  is  the  least  uncertainty  about  where  the  target  will  be  at 
tf.  One  way  to  achieve  a  specified  time  weighting  of  the  guidance  effort  is  to 

take 


tf 

j  N2  dt 

h 

as  the  sub-criterion  to  be  minimized,  where 

x(t)  s  tf  -  t  =  time-to-go  at  t  2:3 


Then: 


if  r  >  0,  resource  expenditures  early  in  the  engagement  are  penalized 
more  than  expenditures  late  in  the  engagement 


if  r  =  0,  there  is  no  time  weighting  of  the  guidance  effort  penalty 

if  r  <  0,  resource  expenditures  late  in  the  engagement  are  penalized 
more  than  expenditures  early  in  the  engagement 

Accordingly,  we  take 

Gost^’N^  2:4 

Ultimately,  means  should  be  sought  to  select  r(tj)  as  a  function  of  the 
estimated  uncertainty  at  tj  regarding  the  present  states  of  the  projectile  and  its 
target  and,  in  particular,  the  predicted  future  states  of  the  target.  Further,  r(tj) 

should  be  chosen,  if  possible,  with  rigorous  consideration  of  the  expected 
statistical  moments  of  errors  in  these  predictions. 

Projectile  and  target  state  estimation  and  the  identification  of  target 
maneuver  strategies  are  not  treated  here,  but  these  become  very  important 
topics  as  the  work  proceeds. 

We  do  not  now  treat  the  case  in  which  r  changes  with  t  after  time  tj, 
because  doing  so  significantly  complicates  the  derivation. 

2.2  Guidance  Error  Penalty 

Consider  a  Cartesian  inertial  coordinate  system  in  which  the  positions 
of  objects  are  described  by  components  x^,  X2,  and  X3.  (See  further  discussion 

in  the  Appendix.)  Then  we  may  write  a  second-order  Taylor  series  to  express 
the  predicted  final  position  of  the  target  for  inertial  axis  k  (k  =  1,  2,  3): 


where  ( )  denotes  an  estimated  or  predicted  value  of  a  variable  and  Ax-p  is  the 
effective  (but  not  a  priori  known)  increment  to  target  acceleration  that,  treated 

A 

as  shown  in  Eq.  2:5,  will  cause  Xpj^(tf)  to  be  an  accurate  estimate  over  the 
interval  of  prediction.  Estimation  of  Ax-p,  an  important  topic,  will  not  be 
dealt  with  in  this  report.  For  the  present,  we  take 

XT  s  xj  +  Axj  2:6 

whence 

A  X  ^  ' 

XTj^(tf)  =  Xp^(t)  +  X  XTj^(t)  +  —  Xpj^ft)  2:7 

For  the  projectile 

X  ^ 

Xj^(tf)  =  Xj^(t)  +  X  Xj^(t)  +  -J-  Xj^(t)  2:8 

provided  the  neglected  higher-order  terms  sum  to  zero  as  t  approaches  tf.  In 
general,  for  large  x,  the  proviso  is  violated.  Nevertheless,  it  is  admissible  to 
consider  a  simplistic  predictive  error  function 

ek(t)  s  ~  ~ 

which  says,  given  x  and  measures  of  the  target  and  projectile  states,  Xp,  Xp,  x,  x; 

given  an  estimate  of  the  effective  target  acceleration  over  the  prediction 
interval  (t,  tf);  and  given  a  particular  value  of  projectile  acceleration  over  the 

same  interval,  a  simplistic  prediction  of  miss  distance  along  axis  k  is  as  stated 
(Eq.  2:9).  We  note  that  the  simplistic  predictive  error  function,  e(t),  is  not  the 
same  as  the  expected  final  miss  distance,  which  we  intend  to  make  essentially 
zero  via  a  two-point  boundary-value  (TPBV)  solution.  Rather,  e(t)  is  used  to 
communicate  to  the  optimization  a  penalty  for  not  nulling  the  simplistic 
prediction  at  every  time  t  (tj  <  t  <  tf). 
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There  is  an  important  difference  between  (a)  use  of  the  simplistic 
predictive  error  function  as  a  sub-criterion  in  a  variational  TPBV  solution 
and  (b)  conventional  guidance  methods  that  act  to  null  some  type  of 
simplistic  prediction  throughout  the  engagement.  The  proposed  optimum 
TPBV  solution  achieves  specified  end  conditions  and  produces  a  trade-off 
between  resource  consumption  and  nulling  of  the  simplistic  prediction  so  as 
to  achieve  zero  terminal  error  and  minimize  the  expenditure  of  resources  in 
the  presence  of  uncertain  information.  In  particular,  this  may  be  done  for 
given  levels  of  uncertainty  about  the  target  maneuver  strategy.  As  limiting 
cases,  for  very  large  uncertainty,  Axj  may  be  taken  as  zero,  and  for  very  little 
uncertainty,  Axj  may  be  assumed  to  be  known  perfectly. 

The  predictive  error  penalty  is  introduced  with  the  function 
3 

k  =  l 

in  which  s(tp  has  a  role  analogous  to  that  of  r,  except  that  s  >  0  corresponds  to 
emphasis  on  the  error  penalty  late  in  the  engagement,  and  s  <  0  pertains  to 
emphasis  early  in  the  engagement. 

2.3  Penalty  on  Loss  of  Kinetic  Enerev 


The  kinetic  energy  of  the  projectile,  per  unit  of  mass,  is 
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k  =  l 


The  time  rate  of  change  of  kinetic  energy  is,  therefore 


X  ^k^k 

k  =  l 
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Thus 


f  ^ 

-  J  Z  ^k^kdt 

ti  k  =  l 

represents  the  loss  in  kinetic  energy  over  the  (tj,  tf)  interval.  The  greater  this 
loss,  the  greater  the  sacrifice  in  range  and  maneuverability  of  the  projectile. 
To  penalize  kinetic  energy  loss,  one  may  introduce  integrals  of  the  following 
sub-criteria 


Gi  =  x-i, 


G2  =  -X2X2, 


G3  =  -X3X3, 


2.4  Limit  on  Total  Number  of  Squibs  Available 

An  approach  toward  limiting  the  total  number  of  squibs  fired  is  to 
introduce  the  penalty 

G4  s  N2  2:12 

2.5  Distance-Travelled  Constraints 

It  may  be  useful  to  maximize  the  distance  travelled  by  the  projectile. 
The  downrange  distance  traversed  is  the  integral  of 

G^  =  sin  Vi  +  X2  cos  Vi  2:13 

and  the  crossrange  distance  is  the  integral  of 

G7  H  X|  cos  Vi  -  ^2  sin  Vi  2:14 

where  Vi  is  the  velocity  heading  at  time  tj.  The  change  of  altitude  (not 
necessarily  to  be  maximized)  is  the  integral  of 


Gg  =  X3 
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The  total  distance  (slant  range)  traversed  is  the  integral  of 

^  (xi  -  t  (x;  -  x;i)x;  ^  (X3  -  X3^)X3 

[(xi  -  x,i)'  +  (X2  -  Xjj)^  +  (xj  -  X3jy]V2 

2.6  Anholonomic  Constraints 

From  Appendix  A  and  first  principles,  one  has  the  following 
anholonomic  constraint  conditions; 

f|  s  mV  +  mg  sin  7  +  qS  Cqq  +  (T*)^  (qSC3  +  CqHi)N  =  0  2:17 

f2  =  mVvcos7  +  T*  ^qSCQH3  -  sincp^N  =  0  2:18 

f3  s  mV7  +  mgcosY+  T*(qSCoH5  -  cosip)n  =  0  2:19 

f^  =  V  -  U|  cos  7  sin  y  -  U2  cos  7  cos  i/  -  U3  sin  7  =  0  2:20 

153  V7  +  sin  7  sin  \|/  +  U2  sin  7  cos  y  -  U3  cos  7=0  2:21 

f^  =  Vy  cos  7  -  cos  v  +  U2  sin  \|/  =  0  2:22 

(7  =  -  Xj  =  0  2:23 

fg  =  V2  -  X2  =  0 

fg  =  V3  -  X3  =  0  2:25 

ho  =  -  Vi  =  0  2:26 

fjl  s  U2  -  V2  =  0  2:27 

f^2  s  U3  -  V3  =  0  2:28 


^13  -  "  Vcos7sin\j/  =  0 


2:29 


fj4  =  V2  -  V  cosy  cos  \|/  =  0 
fl5  =  V3  -  V  sin  7  =  0 

where: 


9 

= 

Ft  +  9i 

2:32 

9 

= 

P 

2:33 

Hi 

= 

1 

2:34 

H3 

-  sin  (p  +  C2  cos  (p 

2:35 

H3 

= 

-  P  (Cjsj^  cos  9  +  C2  sin  (p) 

2:36 

H5 

-  CjvJq^  cos  <P  -  C2  sin  (p 

2:37 

H5 

= 

P  (  ^Na  9  -  C2  cos  (p) 

2:38 

q 

s 

ipv2 

2:39 

• 

q 

= 

p  VV 

2:40 

and  the  parameters  assumed  (for  now)  to  be  constant  are: 


mass 

gravity  acceleration 
atmosphere  density 
reference  area 

zero-normal-force  value  of  drag  coefficient 


induced-drag  factor  (see  Appendix  B) 
normal  force  coefficient 


ly  exp  - 


lyyCOi 


^Ci/2  - 


*From  Hutchings,  Thomas  D.,  Ref.  2,  1984,  which  defines  the  u^rminology  in 
this  factor. 
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(2v)  =  normalized  Magnus  force  coefficient 


scaled  thrust  =  T 


rolling  rate 


squib  burn 
decision  interval 


=  T5t/At«T 


roll  attitude  at  t: 


r'S  9m  15^  “T^  SS  cess  9!Sa  £??1 
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3.  DERIVATION  OF  GOVERNING  EQUATIONS 

3.1  Necessary  and  Sufficient  Conditions  of  the  Calculus  of  Variations 

The  calculus  of  variations  necessary  and  sufficient  conditions  for  existence  of 
a  minimizing  function  F(t)  are  as  follows.  For  the  zeroth-order  freedoms  of  the 
solution,  here  represented  by  u: 


Fu  =  0 


(tj  <  t  <  tf) 


Fuu  ^0  (ti  ^  t  <  tf) 


Euler-Lagrange 

Legendre 


and  for  the  first-order  freedoms,  here  denoted  by  x; 


Fx  =  Fx 


Fx(tf)  =  0 


Fide  ^  0 


(tj  ^  t  <  tf) 


(tf  ^  t  £  tf) 


Euler-Lagrange 


Legendre 


At  least  one  partial  derivative  of  the  form  of  Ineqs.  3:2  and  3:5  must  be 
positive,  and  all  cross  partials  at  the  u  and  x  levels  (such  as  Fy|u2  ^uixi^ 

be  zero  (Kirk,  Ref.  4).  Furthermore: 

F(tf)  =  0  3:6 

3.2  Integrand  of  Performance  Criterion 

From  Eqs.  2:2,  2:10,  2:11,  and  2:12  -  2:15  (omitting  2:16  at  this  time),  the 
integrand  of  the  proposed  performance  criterion  is 

3  3 

F  a  Wq  N^  _  Wj^  XY  x'k  +  W4  n2  +  W5  ^  e^^ 

k  =  1  k  =  1 

+  (x|  sin  Vj  +  X2  cos  +  W7  (x|  cos  Vj  -  X2  sin  Vj/j) 


+  Wg  V  sin  Y  +  X 


1  =  1 
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j  W5  T-S+''  U3  =  W3V3  +  W5  .  X3  -f  X  (ir,  -  V3)  f  Y  XT3 

+  X4  sin  Y  +  ^5  cos  y  -  Xj2  3 

Ineq.  3:2  for  the  zeroth-order  freedoms  provides: 


y  ••  I 


J  W5  X-S+4  ^  0 


which  is  satisfied  provided  W5  is  non-negative. 

Eq.  3:3  applied  for  the  Vp  V2,  V3  freedoms  leads  to: 

^10  “  sin  \|/j  -  W7  cos  Vi  -  ^7  - 


?i^3  3:13 


^11  =  W2  U2  +  2W5  T  ^■^^62  -  Wgcosyj  +  W7sinyj  -  X8  -  X,4  3:14 
^12  ^^3  ^3  2W5T~S'^^  63  +  Wg  -  Xg  -  3:15 

where  ej^  (k  =  1,2,  3)  is  given  by  Eq.  2:9. 
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I*'- 


I 

« 


From  Eq.  3:4  applied  for  V^,  V2,  V3: 


^lOf  ~  ~  ^I2f  -  ® 


>? 

1 


Ineq.  3:5  gives; 

^ViVi  =  I^V2V2  =  ^V3V3  =  0 
and  therefore  these  conditions  are  always  satisfied. 
Next,  Eq.  3:3  applied  for  Xp  X2,  X3  yields: 


Xy  =  2W5T-Se^ 


Xg  =  2W5  X  s 


•I* 

S' 


Xg  =  2W5T-Se3 


3:20 


and  from  Eq.  3:4: 


“hf  -  ^ 


Ineq.  3:5  becomes: 


F**  “F**  sF**  —  0 

xiX]  ”  ^X2X2  ^X3X3  ^ 


For  the  V,  \|/,  and  y  freedoms,  Eq.  3:3  produces: 

>.,m  +  ^4  =  pVS  (x,  [Cd„  +  (TO^CjN]  +  +  >.3H5)Co^'} 

+  (x,2m  +  ^cosY  +  (X3m  +  X5)  y 

-  (Xj3  sin  V  +  Xj4  cos  cos  y  -  Xj5  sin  y  3:23 

^2"^  +  ^6  =  ^6)(V'fsiny  -  V  cosy) 

-  (^4  COS  y  -  X5  sin  y)(Ui  cos  \|/  -  U2  sin  V|/) 

+  X^(U|  sin  Y  +  U2  cos  v)  J  -  ^^3  cos  \|/  +  ^^4  sin  v  3:24 

X3m  +  X5  =  -^  [^  -  (X3m  +  X5)  V  +  mg  (x^  cos  y  -  X3  sin  y) 

-  (X2  m  +  X^)  V  Y  sin  y  +  (X4  cos  y  +  X5  sin  y)  U3 

+  (X4  sin  y  +  X5  cos  y)(u^  sin  \|/  +  U2  cos  v)  J 

-  (X|3  sin  V  +  ^14  cos  \)/)  sin  y  -  Xj5  cos  y  3:25 

•  •  • 

Typically,  X4,  X5,  and  Xg  are  taken  as  zero  from  tj  to  tf.  (See  discussion  in  Section  4.) 
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Eqs.  3:23  -  3:25  involve  V,  y,  and  y,  which  are  functions  of  N  (see  Eqs.  2:17  - 
2:19),  and  Up  U2,  U3,  which,  in  this  case,  are  the  "command"  values  of  projectile 

accelerations  governed  by  the  variational  relationships  in  Eqs.  3:9  -  3:11.  Thus: 

=  ^10  ^  ^11  ^  3;26 


^2  —  ^20  ^  ^21  ^ 


^3  “  -^30  -^31  ^ 


in  which: 


V 

n, 


AlO  *  ^  [pvsx,Cdo  - 

-  (X|3  sin  V  +  cos  v)  cos  y  -  ^15  sin  yj 


“  in^  |2qSm[X|T*C3  +  Co(X2H3  +  X3H5)J 

-  (^2  m  +  X^)  (qS  Cq  H3  -  sin  9) 

-  (X3m  +  ^5)(qSCoH5  -  cos  9  ))  3:30 

^20  =  ^  [(^2"^  +  h)  (^5  -  ^4  cos  y)  (u^  cos  V  -  U2  sin  v) 

+  Xg  (U|  sin  \[r  +  U2  cos  \j/)  -  X13  cos  9  +  X^4  sin  9  J  3:31 

^21  =  '^2^~ofy  (^2^  ^  ^6)  [t*  (qS  C3  +  CoH,)cosy 

-  (qS  Cg  H5  -  cos  9)  sin  yj  3:32 
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which  uses  the  definitions: 


A7  s  X7  +  X,^3  +  sin  Vi  +  W7  cos  Vj 


■^8  ~  ^  ^14  ■*■  ^6  Vi  +  W7  sin  Vi 


Ag  =  X9  +  -  Wg 


From  Eq.  3:4: 


X.|£  m  +  —  0 


(X2£n\  +  X6£)v£C0SYf  =  0 

(^£m  +  ^5f)Vf  =0 


whence: 


A.j£  —  —  X,4f  /m 


4f 


^2f  “  ”  ^6f 


^3f  “  “ 


'5f 


Continuing,  for  the  N  freedom 


Fn  =  2W4N 


and 


=  2WoT''N  +  >,,(T*)2(qSC3  +  C(,H,)  +  Xj  T* 

+  X3T*(qSCoH5  -  COS  cp) 


23 


<  ii.«  ffl.»  tv*  •*•*••*  i'*’i'*"i'i*i^**t*»'i**.**i  *!• 


2  Wo  (t^  N  -  r  t N) 

+  X,  (T‘)2(qSC3  +  CoH,)  4-  X,  (T-)2(qSC3  4  CoH,) 

4  X2T‘(qSC„H3  -  sin  <p)  +  ^2  T*(qS  Cg  H3  +  qSCgHs  -  P  cos  cp) 
4  ^T*(qSCo  H5  +  COS(p)  +  X3  T*(qSCoH5  4  qSCoH;  +  P  sin  cp) 


and  Eq.  3:3  yields,  upon  substitution  of  Eqs.  2:17, 2:40,  and  3:26  -  3:28: 
AjN  +  AjN  +  AgN  =  B 


where: 


A,  2  2WnX* 


A^  2  -  2Wgrxr-l 


An  2  -  2W„ 


3:55* 


B  2  -  T*[T*(qSC3  +  CgHi)Aig  +  (qSCgH3  -  sin(p)A2g 
+  (qSCgHg  -  COS(p)A3g 

~  pvs  (XiT*C3  +  X2CgH3  +  X3  CgH5)(gsinY  +  qSCDg/m) 

+  X2P(qSCgH5  -  cos(p)  -  X3P(qSCgH3  -  sin  cp) J  3:57 

•  • 

The  derivative  N  is  zero  except  at  a  finite  number  of  points  of  singularity. 
Eq.  3:53  must  be  satisfied  between  these  singular  points,  i.e.,  for  the  open  interval 
(tj,  tj  +  At)  and  succeeding  decision  intervals.  Thus 

N  =  (B  -  AgN)/A|  3:58 

From  Eqs.  3:4, 3:47  -  3:49,  and  3:51: 

*  Note  that  all  other  terms  in  N  that  appear  in  Eq.  3:52  cancel  one  another. 
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•  I 


^4f  = 


-  2WoX^J’N^m/T-' 


-  (flfSCoHjf  -  sinipf)  -  X5f(qfSC„H5f  -  cos(p,) 
'  T*(qfSC3  +  CjH,,  ) 

From  Ineq.  3:5: 

A,  =  2WnX^  ^  0 


which  is  always  satisfied  if  Wq  >  0. 

It  can  readily  be  shown  that  all  cross  partials  at  the  "u"  and  "x”  levels  are  zero 
for  the  subject  performance  criterion. 

Lastly,  the  necessary  condition  of  Eq.  3:6  provides  (note  that  the  f  j's  are  zero 
in  Eq.  3:8): 


-  X  *  Wstf-s  X 


k=  1 


k  =  1 


^6  (y\{  Vi  +  V2£  cos  Vi)  +  W7  (Vi£  cos  Vi  -  V2£  sin  Vi) 


+  Wg  V3£  =  0 


Thus,  taking  the  desired  range  of  s,  s  >  0,  requires  that  ^  ej^^  approach  zero 

k  =  1 

quickly  than  Tf“®,  namely 

3 

t  tf  °  3 

In  addition,  if  r  >  0 
3 

(ii)  -  ^  WijVkfUk,  +  W,Nf2  +  W4(v,fSmVi  +  VjfCosVi) 


more 


k=  1 


+  W7(V  l£COSVi  -  V2£sin  Vi)  Wg  V3£  =  0 


If  r  <  0,  it  is  necessary  to  satisfy  3:62  and  3:63  as  well  as 
(iii)  Nf  =  0 

In  the  case  where  r  =  0,  it  is  necessary  to  satisfy  3:62  as  well  as 


3:64 


(iv)  WoNf2  -  ^  WkVkfUkf  +  W4Nf2  +  W6(Vi^sinVi  +  cos  Vj) 


k  =  1 


+  W7  (Vj£  COS  Vi  -  V2£  sin  Vi)  +  ^3  V3£  =  0  3:( 


Solution  of  the  governing  equations  is  discussed  in  the  next  section. 
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4.  SOLVING  THE  GOVERNING  EQUATIONS 

4.1  Establishing  the  Final  Velocity 

To  establish  the  final  velocity  for  off-line  reverse  integration  (data  base 
generation),  if  r  >  0,  Eqs.  3:63  and  3:65  require  that 


Vif(w,v,f  -  Ws  sin  Vi  -  W,  cos  Vi) 

+  V2£(W2V2f  -  cos  Vi  +  WjSinVi) 

+  V3,(W3  V3f  -  Wj)  -  W4  N,2  .  0  4a 

Eqs.  3:1  and  3:9  -  3:11  become ,  at  tf. 

Wi  Vif  =  -  (>,4£  cosYf  -  Xgf  sin7f)sin\i/f  cos  Vf  +  >-iof  4:2 

W2  V2£  =  -  (X4£  cosYf  -  X5£  sinyf^cosVf  +Xg£  sin^f  +  4:^ 

W3  =  “  ^4f  sin  Yf  -  X5£  cos  Yf  +  >>u{  4:4 

However,  Xjq£  =  X|j£  =  Xj2£  =  0  (Eq.  3:16),  and  if  Wj,  W2,  W3  are  not  zero, 
Eqs.  4:1  -  4:4  provide 

[(X4£  cos  Yf  -  ^5£  sin  Yf)  sin  Vf  +  X6£  cos  Vf ]  (v,£  -  sin 

-  W71  cos  Vj) 

+  [(X4£  cos  Yf  -  X5£  sin  Yf)  cos  \g£  -  Xg£  sin  V£]  (V2f 

-  Wg2  COS  Vi  +  W72  sin  Vi) 

+  (X4£  sinYf  +  X5£  cosYf)  (V3£  -  W83)  +  W4  Nf^  =  0  4:5 


where: 


Wfi/Wi,  W71 
W8/W3 


=  W7/W„  W62  =  W6/W2,  W72  =  W7/W2, 


I 


I 


Moreover,  the  inverse  forms  of  Eqs.  2:20  -  2:22  are: 

V,  .  u,  =  (v  cos  y  -  V  y  sin  y)  sin  v  +  V  \Jr  cos  y  cos  v 

Vj  =  Uj  .  (v  cos  Y  -  V  Y  sin  y)  cos  \|/  -  V  vj;  cos  y  sin  \|r 

V3  =  U3  =  V  sin  Y  +  V  Y  cos  Y 
If  N  =  0,  Eqs.  A;22  -  A:24  provide: 

V  =  -  gsiny  -  qSCDg/m 

y  =0 


COSY 

V 


and  Eqs.  4:7  -  4:9  become: 

Vj  =  -q(SCDQ/m)  COSY  sin  y  ' 

V2  =  -  q  (SCpp/m)  COSY  cosy  >  (if  N  =  0)  4:13 

V3  =  -  g  -  q(SCDo/m)sinY  > 

whence  Eq.  4:5  yields 

qf  (SCop/m)  (jifCOSYfSinyf  +  J2f  cos  Yf  cos  yf  +  JafSinyf) 

+  Jlf  (^61  sin  yj  +  W71  cos  yj  +  J2f  (w^z  cos  yi  -  W72  sin  yj) 
+  h{  (Wga  +  g)  -  W4  Nf2  «  0  4:14 


where: 


Jlf  =  (^4fCOSYf  -  ^sinYf)smyf  +  X^fCosyf 
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hi  ■  (^4f  cos  Yf  -  Xgf  sin  Yf)  cos  Vf  -  sin  Vf 


hi  *  ^4fSinYf  +  >-5fCOSYf  4:17 

Now,  substituting  Eq.  2:34,  2:35,  and  2:37  into  Eq,  3:59,  one  obtains,  if 
0  and/or  if  r  >  0, 


[x5fqfSCoC2  +  ^6f(^^^0^Na  ^)1 

T*  (qf  SC3  +  Co) 

(‘if  S  Cg  Cjyj^  ■’■  0  "  ^6f‘lf^^0^2] 
■"  T*  (qfSCa  +  Co) 


smcpf 


cos  (Pf 


For  the  case  <pf  =  0,  Eq.  4:18  becomes: 


^4f  * 


in  which: 


Ksf  =  qfSCgCisjQj  +  1 
^6f  “  ‘If  S  Cg  C2 

Kyf  H  T*(qfSC3  +  Cg) 


whence: 


Jlf  =  (^5fSinVf)|^cosYf  -  sinYf 


+  X^f  cosYf  -  j^cos  Yf  sinxj^f 
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h{  =  V,^5f«>sVf;|^j^cosYf  -  sinYf 


-  Xg,  sinVf  +  j^cosYfCOsv/f 


l3f  =  +  “STfJ- 

and  Eq.  4:14  takes  the  form 

sin  Vf)(K5f  cos  Yf  -  sin  Yf) 

+  hi  (^Of  Vf  -  Kgf  cos  Yf  sin  Vf)  j  |^qf(S  Cog/m)  cos  Yf  sin 

+  W51  sin  Vi  +  W7J  cos  Vj  J 
+  [(>5fCOSVf)(K5fCOSYf  -  Ko£SinYf) 

-  hi  (*^f  sii'  Vf  +  ^6i  Yf  cos  Vf)]  [qf(S  C^p/m)  cos  Yf  cos  Vf 

+  Wg2  Vi  -  W72  sin  Vi  J 

+  [^(KsfSinYf  +  KofCosYf) 

-  XgfKgfSinYf]  [qf(SCDp/m)sinYf+  Wg3  +  g] 

-  Kof  W4Nf2  =  0 

Thus,  from  Eq.  4:26,  defining 

Qf  s  qfS 

one  has  (if  Nf  =  0,  if  r  >  0,  and  if  Vf  =  9f  =  0) 
afQf2  +  bfQf  +  Cf  =  0 


where: 


af  =  Co(Ci3Q/m) 


4:29 


=  (^5f^Na  "  ^6f  ^2)^0  (^83  ■*■  8)sinYf  +  ^sfCoo/n^ 

+  T^C3  [^^83  +  g)  cosYf  -  W4Nf2^ 

+  (Wgi  sin  Vi  +  W71  cos  Vi)  [^(Cq  <^Na  -  T*  C3  sin  yf )sin  Vf 

-  X,g£  (Cq  C2  cos  Yf  sin  Vf  -  T*  C3  cos  Vf)  j 
+  (Wgj  cos  Vi  -  W72  sin  Vi)[^f(Co  cos  Yf  -  T*  C3  sin  Yf)cos  Vf 

-  A.^£(Cq  C2  COS  Yf  cos  Vf  +  T*  C3  sin  Vf)  J 


4:30 


Cf  =  ^5f  (^83  +  g)sinYf  +  T»Co  [x5f(W83  +  g)cosYf  -  W4Nf2] 
+  (w^j  sin  Vi  +  W7J  cos  Vi)  j^^f  (cos  Yf  -  T*  Cq  sin  Yf)  sin  Vf 

+  %T*CoCOSVf] 

+  (w^2  cos  Vi  -  W72  sin  Vi)  (cos  Yf  -  T*  Cq  sin  Yf)  cos  Vf 

-  X5£  T*  Cq  sin  Vf  J  4:31 


Thus 


bf  r/  W  \  2  Cf-i  1/2 

Qf=  -IjJ  (if  Nf  =  Vf  =  ,p,  .  0)  4:32 


If  the  minus  sign  is  used  ^  head  of  the  radical  when  X5£,  Wg,  W7,  Wg,  and  T* 
are  zero,  one  obtains  Qf  =  -  (mg  sin  Yf)/CDQ,  which  requires  that  Vf  =  0  (see 


wMjw  wwulw  wwwufV'wiHiRf  wiivrm*vjtv."v^.»v  <•J^f.^r^’JVJV'^r^rJV)^ 


Eq.  4:10),  and  this  is  not  generally  possible.  Therefore,  we  use  the  plus  sign. 
From  Eqs.  4:27  and  4:32  it  follows  immediately  that 


•<  ■ 


4:33 


4.2  Velocity  Transformations 


The  kinematic  relationships  between  velocities  expressed  in  the 
inertial  (xj,  X2,  X3)  axes  and  wind  (stability)  axes  are: 


V  =  |(Vi2  +  V^2  +  V32)1/2| 

y  =  arcsin  (V3/V) 

and  if  V2  ^  0 

(-^) 

W  cos  y/ 


4:34 


4:35 


y  =  arcsm 
whereas  if  V2  <  0 

y  =  n  -  arcsin 


4:36 


in(^) 

VV  cos  7' 


4:36 


If  7  =  ±  71/2  (i.e.,  if  cos  7  =  0),  y  is  undefined.  In  Eq.  4:35,  7  varies  between  ± 
7t / 2,  and  in  Eqs.  4:36,  y  ranges  between  -  7t/2and  +  3n/2. 


4.3  Estimation  of  Time-to-Go  (t) 

A  simplistic  estimate  of  time-to-go  is 
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Closing  Velocity 

AV 

cos  0 

3 

k  =  1 

COS  0 

.V 

where  the  j's  are  unitless  direction  vectors,  0  is  the  angle  between  the  vectors 
Ax  and  AV,  and: 


-  ""k 

^^k  =  Xxk  -  Vk 


cos  0  = 


Ax  AV 


Ax  AV 


''  --.L^ _  4-41 

t  3 

X  Axi^  AVk 
k  =  1 

Further  discussion  of  to  time-to-go  estimation  is  presented  in 
Appendix  D. 


4.4  Firing  Commands  for  Multiple- Sector  Thrusters 

The  CAT  projectile  may  be  designed  with  multiple  (perhaps  four) 
sectors  of  impulsive  thrusters,  each  sector  containing  a  number  of  squibs  that 
fire  through  a  common  sector  orifice.  For  any  given  sector  there  is  an  initial 
roll  attitude,  (pj:,  the  index  j  denoting  the  sector.  Correspondingly,  Eqs.  A:29  - 
A:30  and  2:35  -  2:38  provide  multiple  values  of  H3j,  H5j  as  functions  of 

time.  As  a  result,  Eq.  3:30  provides  a  different  ^|j(t)  value  for  each  sector,  and 
Eq.  3:57  and  3:58  give  a  different  Nj(t)  value  for  each  sector. 

The  other  equations  of  the  OPTG  solution  are  the  same  for  all  sectors. 
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4.5  Numerical  Solution  Procedure 


To  start  numerical  integration  of  the  differential  equations,  the 
procedure  is  as  follows: 

(1)  Describe  projectile  (m,  S,  C^q,  K,  Cjsj^,  and  lesser  constants  Ig,  Ij,  Cp 

coi,  0)2,  lyy,  d^ef,  Vgy,  '^rnuzzle'  environment  (g,  p),  and  target 

(xj.j  (t),  Xj2  (t),  xjg  (t)). 

(2)  Specify  guidance  parameters  (Wq,...,  Wg,  r  >  0,  s  >  0). 

(3)  Specify  P  and  At  (At  =  0.05  sec.).  Compute  firing  constants  T*,  Cg,  C2, 
and  C3  from  Eqs.  A:17,  A;19,  A:31  and  A:32,  respectively. 

For  Forward  Integration 

(4.F)  Set  x-(  =  X2  =  X3  =  0  and  N  =  0.  Specify  initial  \j/,  y,  and  (p.  Set  V  = 
^muzzle-  Specify  Xi,...,Xi5  at  tg. 

For  Reverse  Integration  (Negative  At) 

(4.R.1)  SetX|  =  xti£,  X2  =  xt2£,  X3  =  xt3£.  Specify  N£,  yp  Xgp  ?ig£.  Use\|/£and 
(p£  of  zero. 

(4.R.2)  Compute  V£  from  Eqs.  4:29  -  4:33,  adjusting  W4,  N£,  A.5,  and  if 
necessary  to  obtain  acceptable  V£. 

(4.R.3)  Compute  X4£  from  Eq.  4:18  emd  from  Eqs.  3:47  -  3:49.  Set 

=  Xg£  =  Xg£  =  A.|Q£  =  Xjj£  =  Xj2£  =  0,  per  Eqs.  3:21  and  3:16. 

Now  the  system  states  are  initialized  (or  "finalized"),  and  the  state 
derivatives  may  be  calculated  as  follows: 


•  <  *,1  4,4  ♦  J  1.4  I 


(5)  Compute  Vj,  V2/  V3  from  Eqs.  2:29  -  2:31.  Compute  ej,  e2,  e3  from  Eq. 

2:9.  Compute  range  from  ((xyi  +  (xj2  -X2)^  +  (x-j-3  -X3)2)^/2. 

A  A 

(6)  If  range  =  0,  setx  =  0;  otherwise  compute  t  from  Eq.  4:41. 

A 

(7)  If  X  =0,  compute  Uj,  U2,  U3  (projectile  accelerations)  from  the 
kinematic  relationships,  Eqs.  4:7  -  4:9;  otherwise  compute  Uj,  U2,  U3 
from  the  variational  solution,  Eqs.  3:35  -  3:37. 

A  •  * 

(8)  Ifx  =  0,  setNf  =  0;  otherwise,  if  unused  squibs  remain,  compute  N 
from  Eqs.  3:55  -  3:58. 

(9)  IfN  >  Threshold,  set  N  =  1.  (Fire  squib.) 

If  N  ^  Threshold,  set  N  =  0.  (Don't  fire  squib.) 

•  •  • 

(10)  Compute  wind  (stability)  axes  rates  V,  \|/,  y  from  Eqs.  A;33  -  A:35. 

(11)  Compute  Xy  ^2/  ^3  from  Eqs.  3:26  -  3:34. 


•  •  « 


(12)  If  X  =  0,  set  =  X3£  =  X^^  =  0;  otherwise  compute  Xy,  Xg,  Xp  from 
Eqs.  3:18  -  3:20. 

•  •  • 

(13)  Compute  ^12'  using  Eqs.  3:38  -  3:40. 

This  completes  the  calculations  of  state  derivatives  from  state  values. 
Numerical  integration  may  now  be  employed  indefinitely  to  compute  the 

states  ((p,  Xj,  X2,  X3,  V,  xp,  y,  and  Zi)  from  their  initial  (or  final)  values  and  their 

•  •  •  • 

respective  derivatives  (P,  V|,  Vy,  V3,  V,  \|r,  y,  and  U. 

4.6  Discussion 

The  best  decision  interval  (At)  is  a  function  of  projectile  rolling  rate  (P). 
Preliminary  evidence  suggests  that  the  difference  P  -  1/At  should  be 
approximately  +  1.0  Hz.  This  has  the  effect  of  providing  roll  attitude  (9) 


o 

i 

i 

i 

vM* 


w 

o 

I 

n 

•iw 

i 

i 

o 


I 


segments  of  approximately  18  deg.  between  20  discrete  angles  "available"  for 
thruster  firing,  with  each  angle  becoming  available  once  per  second. 

The  N  command  equation  (Eq.  3:58)  may  be  used  in  formal  integration 
of  the  equations  of  motion  and  adjoint  equations,  proceeding  from  specified 
initial  physical  and  adjoint  (variable  Lagrange  multiplier)  states,  or  it  may  be 
used  in  reverse  integration  from  specified  final  (burst  or  impact)  conditions. 
The  use  of  reverse  integration  is  helpful  in  mapping  the  behavior  of  the 
system  as  a  function  of  the  constant  Lagrange  multipliers  (the  W's),  ^5,  and 
Xg.  When  employing  reverse  integration,  final  values  for  the  adjoint  states 
are  known  from  the  variational  conditions,  except  that  X5£  and  values 

must  be  stipulated  from  consideration  of  the  required  initial  physical  states. 
(Here  we  assume  ^4  =  ^5  =  ^6  ~  choose  to  refer  to  X5,  as  the 

"steering"  multipliers,  because  they  are  the  principal  mechanisms  for 
satisfying  the  two-point  boundary  values  imposed  on  guidance  of  the 
projectile.  (The  other  ^’s  and  the  W's  also  influence  the  boundary  values,  but 

serve  primarily  to  determine  the  optimal  character  of  the  solution. 

» 

Furthermore,  the  final  values  of  most  of  the  other  ?t’s  and  the  X's  are  known 
immediately  once  and  X(^^  have  been  provided. 

«  *  • 

In  reverse  integration,  ej,  62,  63  and  Xy,  X^,  Xg  begin  at  zero.  In  view  of 

•  • 

Eqs.  3:21,  X7,  Xg,  X9  also  begin  at  zero.  As  long  as  e^  07, 63  remain  zero, 

X-^2  (ignoring  W.j,  W2,  W3,  W^,  W7,  and  Wg)  will  be  proportional  to  X14, 

•  •  • 

respectively.  Now,  if  X,jq,  Xjj,  are  non-zero,  A,jq,  X|2  will  integrate 
away  from  their  beginning  zero  values  at  tf  (Eqs.  3:16).  As  this  occurs.  Up  U2, 
U3  will  become  biased,  which  is  inappropriate  to  the  purpose  of  maintaining 
ep  e2,  eg  at  zero.  We  conclude,  therefore,  that  it  may  be  desirable  that  X^g,  ^^4, 
X|5  be  kept  equal  to  zero  at  all  times. 

The  W5  G5  constraint,  involving  a  simplistic  prediction  of  final  error, 
introduces  the  X7,  Xg,  X9,  X|q,  Xjp  Xj2  adjoint  differential  equations,  coupling 
them  to  the  Xp  X2,  Xg  adjoint  differential  equations  via  the  variational  (i.e., 
the  "command"  as  contrasted  with  the  kinematic)  expressions  for  Up  U2,  Ug. 
As  Wg  0,  the  coupling  becomes  weaker,  disappearing  altogther  at  Wg  =  0. 


Lfki'a.t'l.tv 


If  W5  =  0,  the  ^4,  X,5,  X.g,  Xj|,  and  X-^2  differential  equations  may,  in 
principle,  be  eliminated  from  the  solution,  with  the  present  roles  of  X4,  X5, 
and  then  assumed  by  ^15'  ^14- 


Two  configurations  of  the  operational  OPTG  guidance  should  be 
considered.  In  the  first,  polynomial  networks  would  be  used  to  estimate  the 
initial  values  of  all  the  k's  and  to  re-initialize  the  solution  during  the  flight  of 
the  projectile.  For  example,  before  the  gun  is  fired,  the  fire  control  system 
(FCS)  could  use  polynomial  networks  to  compute  initial  values  of  X^,  X2,  Xy 
with  the  initial  values  of  Xy  ...,  >.^5  being  the  same  for  all  trajectories  and 
therefore  not  requiring  pre-fire  estimation.  Then,  with  the  projectile  in 
flight,  other  polynomial  networks  in  the  FCS  could  be  used  to  modify  X5  and 
X^  to  correct  for  target  maneuvers  and  atmospheric  disturbances  of  the 
projectile  trajectory.  In  the  second  formulation,  polynomial  networks  would 
provide  the  N  decision  (squib  fire/no-fire  decision).  These  configurations 
will  be  described  further  in  future  reports. 

The  OPTG  formulation,  implemented  via  either  of  the  above 
configurations,  will  aid  the  system  as  it  deals  with  information  uncertainties 
and  information  denial  and  will  provide  efficient  use  of  the  projectile 
impulse  thrusters.  Regarding  the  problem  of  uncertainties,  the  OPTG 
guidance,  implemented  using  APNs,  will  (a)  respond  optimally  to  predictable 
maneuvers  of  the  target  and  (b)  optimally  desensitize  the  guidance  to  the 
adverse  effects  of  unpredictable  target  motions  and  tracking  errors.  These 
benefits  will  accrue  because,  with  careful  simulation  and  data  base  generation 
for  APN  synthesis,  the  true  non-parametric  statistics  of  the  problem  will  be 
mastered.  (Incidentally,  the  payoff  for  good  simulations  will  be  measurably 
increased.)  Regarding  the  efficient  management  of  projectile  thrust  impulse, 
the  calculus  of  variations  will  provide  a  data  base  of  variational  extremal 
trajectories  that  optimize  the  usage  of  projectile  impulse,  correcting  by  the 
most  appropriate  numbers  of  squib  firings  as  the  flights  progress. 
Conventional  guidance  laws  do  not  have  the  capability  (the  intelligence)  to 
do  this. 
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5.  A  BASELINE  GUIDANCE  LAW 

If  the  simplistic  predictive  error  function  of  Eq.  2:9  is  set  to  zero,  one 
obtains 


Uck  =  f  [’‘Tk  -  Xk  +  ^  (xtlc  -  4)  +  T’>‘Tk] 


in  which  the  subscript  "c"  denotes  a  command  value.  Thus,  from  Eqs.  2:16  - 
2:18: 

Wf,  =  cos  Y  sin  xff  +  U(;2  cos  y  cos  y  +  sin  y  5:2 

(Vylj;  =  -  Ugj  sin  y  sin  y  -  sin  y  cos  v  +  cos  y  5:3 

(V  v  cos  y)^  =  cos  v  -  Uj,2  sin  >j/  5:4 


Eqs.  A:33  -  A;35  now  become: 


K*(v 


c  +  g  sin  y  + 


M  (1)  =  — : - 

c  (T*)(qSC3  +  CqHO 


qs  Cdq) 


K*  (v  y  cos  y)^ 


N  (2)  =  - 1—1 - lli^ 

*■  qS  Cq  H3  -  sin  9 


Nc^3)  = 


K*  [(Vy)^.  +  g  cos  y] 
qS  Cq  Hj  -  cos  9 


where 


K*  H  -  m/T*  5:8 

In  practice,  K*  may  take  on  a  different  value  from  that  above  so  as  to  provide 
approriate  closed-loop  responsiveness  and  stability  of  the  guidance  system. 
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Usually,  the  first  command  equation,  i.e.,  on  would  not  be 

implemented,  because  control  over  V  is  weak  when  the  directions  the 
thrusters  fire  are  always  normal  to  the  projectile  x  axis. 

t 

The  other  two  values  of  in  Eqs.  5:5  -  5:6  may  each  be  submitted  to 
testing  against  respective  thresholds,  and  the  decision  to  fire  an  impulsive 
thruster  may  be  made  if  one  or  more  of  these  threshold  tests  is /are 
affirmative.  However,  such  a  policy  leads  to  firing  if  the  expected  result  along 
any  axis  is  favorable  in  spite  of  any  adverse  repercussions  along  the  other 
axes;  therefore  a  decision  based  on  a  weighted  sum  of  the  individual  sub¬ 
decisions  is  employed,  taking 

Nc  =  |sin  (p|  +  jcos  <pj  5:9 

The  threshold  on  is  nominally  set  at  zero,  but  may  be  raised  to  account  for 
information  uncertainties. 


Combining  Eqs.  2:28, 5:4,  and  5:6 


K*  (Uij,  cos  V  -  U2^  sin  y) 
qS  Cq  C2  cos  9  -  (qSCoCfj  +  l)sin9 


and  combining  Eqs.  2:30,  5:3,  and  5:7 


K*  [U|^  sin  Y  sin  y  +  U2j,  sin  y  cos  y  -  +  g)  cos  yj 

(qSC(,  ^Na  qSCoC2sin9 


Eqs.  5:1,  5:9,  5:10  and  5:11  define  a  baseline  guidance  law,  diagrammed 
in  Figure  2,  that  may  be  studied  in  comparison  with  the  OPTG  formulation. 
Note  the  equivalence  of  Eqs.  3:9  -  3:11  and  5:1  as  W5  tends  to  infinity.  In  the 

limit,  the  OPTG  guidance  law  responds  in  the  same  way  as  the  baseline  law  if 
total  emphasis  is  placed  on  the  Wg  term  in  the  integral  performance  criterion. 


Figure  2:  Baseline  Guidance  Computations 


There  are  important  reasons  for  not  reducing  the  OPTG  guidance  to  the 
baseline  form.  First,  the  OPTG  guidance  is  a  two-point  boundary-value 
guidance  law;  it  provides  essentially  zero  final  error  if  the  maneuvers  of  the 
target  are  within  allowable  limits  for  the  projectile  system.  On  the  other 
hand,  the  baseline  guidance  law  acts  to  minimize  the  integral  of  squared 
predicted  error;  with  the  baseline  formulation,  the  actual  final  error  receives 
no  more  attention  mathematically  than  does  the  final  error  predicted  at  each 
of  the  decision  points  prior  to  tf.  Ultimately,  it  is  not  the  time  history  of 
predicted  final  errors  but  the  actual  error  at  tf  that  determines  the  probability 

of  target  kill.  OPTG  guidance  addresses  directly  the  requirements  on  final 
error  and,  depending  on  the  value  selected  for  Wg,  may  put  little  or  almost 

no  weight  on  the  time  history  of  predicted  final  errors. 

A  second  important  reason  for  using  the  OPTG  formulation  in 
preference  to  the  baseline  guidance  law  is  that  OPTG  guidance  provides 
conservation  of  projectile  maneuver  resources  in  a  way  that  does  not 
compromise  final  error.  The  Wg,  Wj,  W2,  W3,  W4,  Wg,  W7,  and  W3  terms  in 

the  OPTG  formulation  address  the  management  of  maneuver  resources. 
While  it  is  possible  to  incorporate  a  resource  expenditure  penalty  in  the 
baseline  type  of  guidance  law,  the  effect  of  this  penalty  would  be  to  engender 
greater  miss  distances.  In  other  words,  maneuver  resource  management  is 
best  treated  in  the  two-point  boundary-value  guidance  context. 


|*|.a<4^B*|  a'*  •*•  •«!  a*i'*a«4'a*i'k*A'Jl'A'i 


la* 


4 


<  z' 


6.  SIMULATION  PROGRAM 


The  essential  parts  of  the  three-degree-of-freedom  (3DOF)  simulation 
employed  are  listed  in  Appendix  E.  A  goal  of  the  program  is  to  enable  one  to 
examine  flights  starting  from  arbitrary  points  in  space  and  proceeding  either 
forward  or  backward  in  time  to  suitable  stopping  points.  A  system  state 
structure  is  therefore  defined  with  which  one  can  easily  store  and  return  the 
required  system  "snapshots",  or  sets  of  variables  which  uniquely  describe  the 
status  of  the  engagement  under  study.  The  states  S  =  (tp,  ^  V,  y,  \f,  2iu  xj'  ^T' 

•  «  I  • 

N)  are  arrived  at  by  integrating  their  respective  derivatives,  S  =  (P,  V,  V,  y,  v, 

•  •  • 

L  Yj  U-p,  N)  over  time  (where  N  is  treated  as  either  0  or  1,  as  discussed  in 
Section  4.5). 


The  state  derivatives,  S,  can  be  written  entirely  as  functions  of  S;  in  the 
program,  these  calculations  are  split  into  common  {dyn)  and  squib-related 
(deriv)  parts.  Note  that  the  kinematic  equations  in  dyn  and  deriv  are  general, 
the  variational  equations  derive  from  the  calculus  presented  in  Section  3,  and 

the  required  system-specific  constants  are  isolated  in  an  include  file 
isystem.h).  The  guidance  command,  N^.,  subjected  to  a  user-defined 

threshold  for  the  fire/no-fire  decision,  is  calculated  according  to  one  of  three 
guidance  modes: 


Ballistic  (N  =  0) 

Baseline  (Section  5) 

Optimum  Path-to-Go  (OPTG;  Section  4). 


The  fly  routine  integrates  from  Sj^  to  Sj^^.  j  (or  vice-versa  given  the 

derivatives  at  the  appropriate  starting  point)  through  either  the  rectangular 
method  or  Heun's  method  of  integration,  as  selected  by  the  user.  In  addition 
to  performing  a  single  iteration  of  the  system,  fly  examines  whether  or  not 
the  timestep  should  be  reduced  or  reversed  to  "zero-in-on"  the  desired  final 
condition  (such  as  minimum  time-to-go,  x,  for  closest  point  of  approach). 
Not  shown  in  the  listing  is  the  capability  of  storing  the  system  state  and 
performing  a  forward  look-ahead  —  a  rapid  extrapolation  of  current 


information  to  examine  its  closed-loop  implications  —  which  is  required  if 
miss-distance  sensitivity  APNs  are  used. 

The  system  is  flown  as  long  as  allowed.  The  conditions  allowing  flight 
vary  with  flight  direction  according  to  the  routine  conds_alloxv.  Also 
depending  on  direction,  the  system  is  initialized  or  finalized  in  main  after  the 
state  and  control  variables  have  been  retrieved  from  the  user  interface  and 
translated  into  internal  units.  The  interface  (not  listed)  is  a  user-friendly 
general  package  developed  by  Barron  Associates  to  improve  flexibility  and 
portability  in  simulation  development.  With  it,  the  user  may  save  or  restore 
named  sets  of  variables,  change  as  many  or  as  few  parameters  as  desired 
between  trials,  automatically  reset  values  to  defaults,  and  get  descriptive 
information  on  the  purpose  of  each  program  option. 

When  the  user  quits  a  series  of  simulated  flights,  or  trials,  in  one 
direction  (more_trials),  the  option  to  examine  flights  in  the  other  direction  is 
presented.  In  addition  to  end-point  condition  manipulation  through  the 
interface,  the  user  may  define  a  set  of  parameters  to  adjust,  and  initiate, 
certain  automatic  explorations  according  to  the  following  search  modes  found 
in  advise: 

0:  None:  use  interface 

1:  Grid:  search  in  a  regular  pattern  from  loops  within  loops. 

2:  Grope:  employ  the  Guided  Random  Optimizer  of  Performance 
Error  (GROPE),  a  search  algorithm  developed  by  Barron  Associates 
(not  listed).  GROPE  attempts  to  minimize  a  criterion  defined  by  the 
user  in  main,  (for  ex:  miss  distance).  The  search  algorithm  is 
capable  of  efficient  exploration  of  high-dimensioned  multi-modal 
surfaces. 

3:  Gauss:  randomly  select  values  for  parameters  according  to  a  zero- 
mean  gaussian  distribution  of  user-defined  standard  deviation. 

4:  Uniform:  employ  a  uniform  distribution  bounded  in  each 
parameter  by  the  user. 

5:  "Firtree":  use  a  special  sensitivity-APN  database  generation  mode. 
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The  bounds  to  employ  when  changing  the  parameters  are  defined  in  the 
include  file  bounds.h. 

Also  present  in  the  program  are  options  to  select  a  wide  or  narrow 
screen  of  published  information  separated  by  adjustable  time  intervals, 
output  a  history  of  state  variables  for  presentation  or  analysis,  enter  a  "debug" 
mode  to  receive  interior  details  on  certain  calculations,  define  the  record 
modes  to  employ  for  efficient  database  generation,  and  to  "bovmce",  or  return 
in  the  opposite  direction,  according  to  the  same  or  a  different  guidance  mode. 
(For  instance,  one  could  simulate  a  reverse  trajectory  under  OPTG  guidance 
and  "bounce"  forward  ballistically  to  compare  the  divert,  or  difference  in  final 
position.)  Not  shown  is  the  use  of  adaptively-synthesized  polynomial 
networks  (APNs)  for  guidance  initialization  and/or  adjustment  in  the  closed- 
loop  OPTG  system. 


7.  OPTG  AND  BASELINE  GUIDANCE  SIMULATION  RESULTS 

Use  of  the  3DOF  CAT  projectile  simulation  software  presented  in 
Section  6  shows  the  OPTG  guidance  law  consistently  producing  smaller 
terminal  miss  distances,  using  fewer  squib  firings,  and  providing  greater 
maximum  divert  capability  than  does  the  baseline  guidance.  Section  1 
provides  some  discussion  of  these  effects;  here  we  present  detailed  results 
supporting  the  conclusions  reached. 

7.1  Engagements  Simulated 

The  engagement  geometry  simulated  is  that  of  a  maneuvering  or  non¬ 
maneuvering  target  located  in  the  X|-X3  plane  at  a  distance  X2  =  10,000  ft.  from 
the  gun  muzzle  (where  x^,  X2,  X3  correspond  to  "east",  "north",  height 
position).  Target  accelerations  are  tracked  with  some  degree  of  uncertainty 
(gaussianly-distributed  random  measurement  noise)  and  are,  in  the 
simulations  studied,  either  constant  or  governed  by  the  equations: 

uTi(0  =  g  sin(cot)/3  ft./sec.^  7:1 

073(0  - -32.8  0)2  sin((ot)  ft./sec.2  7:2 

where  O)  is  O.Iti  rad. /sec.  The  equations  are  derived  from  the  elliptic  target 
motion  relationships  defined  in  Reference  4: 

y(t)  =  g  ( t  -  sin(o)t  +  (py)/to  )/3co  meters  7:3 

z(t)  =  50  +  10  sin(o)t  +  9z)  meters  7:4 

where  (y,  z)  is  (x|,  X3),  and  the  random  target  phase  components,  (py  and  cpz/ 
are  omitted. 

Initial  projectile  conditions  were  constant:  angle  of  inclination,  y  =  2°; 
muzzle  velocity,  Vm  =  3821.5  ft./sec.;  roll  attitude,  9  =  0;  muzzle  attitude  east 
of  north,  \y  =  0;  (constant)  roll  rate,  P  =  21  Hz.;  and  decision  time  step,  At  =  0.05 
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sec.  (Note  that  the  latter  two  parameters  cause  an  effective  decision  roll  rate 
of  1  Hz.)  An  unguided  projectile  (i.e.,  ballistic  trajectory)  was  found  to  take 
approximately  3.065  seconds,  and  to  cross  the  10,000  ft.  downrange  plane  at  an 
X|  ("east")  position  of  0.0  ft.,  and  an  X3  (height)  position  of  182.23  ft.  --  a 

location  henceforth  referred  to  as  the  "ballistic  point". 

7.2  Simplified  OPTG  System 

The  complete  OPTG  guidance  system  envisioned  (diagrammed  at  a 
high  level  in  Figure  l.a)  is  detailed  in  Figure  3. a.  Adjoint  differential 
equations  are  integrated,  and  the  squib  fire  command  calculated,  at  each 
decision  point  to  determine  whether  a  squib  will  fire,  diverting  or  not 
diverting  the  projectile  at  that  point  in  its  trajectory.  Projectile  and  target  state 
information  is  filtered  (potentially  using  linear  polynomial  filters,  a  class  of 
estimation  and  prediction  filters)  to  track  progress  and  update  the  OPTG 
guidance  equations.  Those  equations  are  initialized  and,  if  necessary, 
periodically  re-initialized,  by  an  APN  synthesized  off-line  from  simulations 
of  representative  trajectories.  The  APN  of  Figure  3.a  relies  on  a  preceding 
explicit  prediction  of  the  intercept  point;  but,  as  shown  earlier  in  Figure  l.b, 
this  task  may  be  subsumed  by  a  more  comprehensive  APN  which  works 
directly  from  the  filtered  target  and  projectile  state  measurements. 

The  simplified  OPTG  guidance  system  prototyped  and  evaluated  in  this 
first  phase  of  work  is  shown  in  Figure  3.b.  As  the  class  of  engagements 
studied  terminates  at  a  nearly  constant  range,  intercept  point  prediction  is 
based  on  the  ballistic  flight  time.  The  target  position  estimator,  using  this  t 
and  the  true  initial  target  position,  was  simulated  to  be  only  a  minor  source  of 
noise  in  the  experiments.  (A  complete  examination  of  this  important  topic  -- 
beyond  the  scope  of  this  study  -  will  be  required  for  full  Optimum  Path-To- 
Go  guidance  implementation).  Note  that,  unlike  full  OPTG  guidance,  the 
simplified  version  performs  no  guidance  re-initialization  or  adjustment;  it  is 
the  "open  loop"  component  of  the  eventual  system. 

In  the  cases  examined,  the  constant  OPTG  weights  defining  the  basic 
character  of  the  trajectories,  Wg  through  Wg,  were,  in  order:  { 10,  0,  0,  0, 100, 

0.001,  0,  0,  0  ).  Note  that  the  adjoint  equation  weight  which  employs  target 


state  feedback  (W5)  was  disabled;  in  fact,  only  the  Nc  scaling  weight,  Wq,  and 
the  squib  conservation  weight,  W4,  were  used  to  govern  the  basic  OPTG 
tradeoffs.  (The  feedback  weight,  W5,  though  shown  with  a  small  value,  is 
effectively  disabled  when  the  kinetic  energy  weights,  W2,W2,W3,  are  absent.) 
Although  most  of  the  weights  (such  as  W5,  W7,  Wg,  the  downrange, 
crossrange,  and  altitude  change  weights)  were  studied,  it  was  found  that  the 
small  subset  above  was  sufficient  to  demonstrate  the  necessary  concepts. 
However,  further  research  could  likely  uncover  more  trajectory 
enhancements  by  examining  the  full  use  of  these  temporarily  disabled 
weights. 

7.3  OPTG  vs.  Baseline  Guidance 

To  examine  the  operational  range,  or  "handprint",  of  the  optimum- 
path-to-go  strategy,  initial  values  of  two  of  the  controlling  Lagrange 
multipliers,  X2  and  X3,  were  varied  randomly  in  simulation  between  values 

of  +/-  1.50,  with  all  the  other  ^'s  starting  at  0.0.  (Note  that  OPTG  guidance 
reduces  to  the  ballistic  case  when  the  entire  vector,  is  zeroed.)  The  database 
of  engagements  thus  generated  was  filtered  to  reduce  redundant  impacts, 
retaining  trajectories  firing  the  fewest  squibs  in  the  case  of  a  conflict.  The 
Algorithm  for  Synthesis  of  Polynomial  Networks  (ASPN)  was  then  employed 
to  synthesize  a  model  (an  APN)  capable  of  revealing  the  appropriate  ^ 
values  to  use  when  given  the  expected  final  target  position  (xjif  and  xT3£). 

Using  "E"  and  "H"  for  expected  position.  Figure  3.c  details  the  resulting 
initializing  APN,  which  network  consists  of  a  pair  of  "carved  double" 
elements,  surrounded  by  "normalization"  and  "unitization"  nodes  (Ref.  12). 

When  installed  in  the  simulation,  and  relied  on  to  initialize  the  OPTG 
guidance  for  a  host  of  target  position  cases  within  the  maneuver  envelope, 
the  APN  proved  to  be  very  successful.  Against  motionless  targets,  the  OPTG 
handprint  in  the  miss  distance  plane,  shown  in  Figure  4.a,  is  quite  circular, 
providing  good  controllability  over  the  entire  region  physically  reachable  by 
the  projectile.  Examination  of  the  data  reveals,  importantly,  that  within  the 
performance  boundaries,  the  accuracy  of  the  OPTG  system  does  not  degrade 
substantially  with  required  divert,  or  final  distance  from  the  ballistic  impact 


Approximate  Equations: 


lambda_2  =  32.10  -  0.5014*H  +  0.3717e-2*E  + 
0.2668e2'»H''2  -  0.4849e-3*H'^3 

lambda_3  =  -4.919  +  0.07874*H  +  0.01657*E  - 

0.4292e-3*H^2  -  0.7073e-4*E^2  +  0.2963e-5*E^3 
+  0.7799e-6*H^3 


Figure  3.c:  Prototype  Initializing  APN 


Divert  Distance  from  Baiiistic  impact  (ft.) 


Figure  4.b:  OPTG  Guidance  Squib  Usage  (0  Threshold) 


point.  Figure  4.b  reveals  that  there  is  a  clear  relationship  between  squibs  fired 
and  the  resulting  divert;  in  fact,  the  number  of  squibs  needed  in  a  particular 
engagement  is  about  one  half  the  total  divert  (in  feet)  required. 

When  the  baseline  guidance  law  is  sent  after  the  same  set  of  targets,  the 
handprint  depicted  in  Figure  5.a  results.  Its  pattern  is  more  that  of  a  diamond 
than  a  circle;  significantly,  points  at  the  corners  of  the  OPTG  maneuver 
boundary  appear  to  be  unreachable  with  the  baseline  method.  Many  more 
squibs  are  being  fired,  as  shown  in  Figure  5.b,  but  no  clear  relationship  exists 
between  thruster  energy  expended  and  the  resulting  divert  (a  consequence 
partly  of  the  imbridled  continuous  correcting  performed  by  the  baseline  law). 

Note  that  if  the  projectiles  in  the  above  baseline  case  were  limited  to  40 
squibs,  most  trajectories  would  have  had  many  fewer  opportunities  to  refine 
final  accuracy,  vdth  a  consequent  (but  slight)  worsening  of  miss  distance.  The 
difference  is  only  slight  because,  for  the  squibs  to  divert  the  projectile 
significantly,  they  must  be  fired  early  enough  for  the  effects  to  integrate 
throughout  the  flight.  The  loss  of  final  refinements  when  resources  are 
exhausted  is  therefore  a  consideration  less  important  than  the  proper  early 
use  of  those  resources.  Both  issues  are  successfully  addressed  by  the  OPTG 
guidance  method,  which  inherently  conserves  resources  by  acting  early  -  a 
property  not  shared  by  late-acting  pursuit,  proportional  navigation,  or 
predictive  proportional  navigation  (e.g.,  the  baseline)  guidance  laws. 

Putting  a  threshold  on  the  guidance  command  equation  -  essential 
when  uncertainties  in  target  accelerations  exist  --  reduces  the  squib  firing 
called  for  by  the  baseline  guidance  law.  Figures  6.a,b  depict  the  results  when  a 
threshold  of  0.2  is  employed  in  the  baseline  law.  Now  a  definite  relationship 
exists  between  squibs  fired  and  divert  achieved  (6.b),  and  though  the  system 
uses  fewer  squibs  than  without  a  threshold,  it  still  fires  roughly  8  more  squibs 
than  with  OPTG  guidance.  The  price  of  saving  squibs  through  a  threshold, 
however,  is  revealed  in  the  resulting  miss  distance,  as  suggested  by  the 
distorted  and  shrunken  handprint  (6.a).  The  accuracy  of  the  conventional 
guidance  law  in  this,  the  simplest  of  engagement  scenarios,  is  good  below  a 
divert  of  about  32  feet,  as  revealed  in  Figure  7.a,  but  worsens  quite  sharply 


Final  East  Position,  xl  (ft.) 

Figure  6.a:  Baseline  Guidance  Handprint  (0.2  Threshold) 


beyond  that.  When  the  thresholds  necessary  for  use  in  noisy  environments 
are  employed,  the  accuracy  of  the  baseline  guidance  only  worsens  (Figure  7.b). 


The  OPTG  guidance  method  is  only  slightly  affected  by  use  of  the  same 
guidance  command  equation  threshold  (0.2).  Though  the  APNs  were  not 
trained  on  such  a  threshold,  its  presence  served  mostly  to  delay  the  onset  of 
each  fire  decision,  resulting  in  slight  "twists”  of  the  intended  trajectories.  The 
end  effect  was  minor,  and  the  trajectories  have  a  very  similar  miss  distance 
distribution,  confirming  that  OPTG  laws  can  be  made  much  less  sensitive  to 
the  values  of  any  thresholds  employed  to  deal  with  uncertain  (noisy) 
information.  Figure  8  reveals  that,  with  no  squib  fire  decision  threshold,  the 
OPTG  trajectories  achieved  given  divert  distances  with  about  10  fewer  squibs 
than  the  thresholded  (i.e.,  most  competitive)  baseline  method. 

Note  that,  as  initialization  and  re-initialization  of  an  optimum-path- 
to-go  guidance  law  occurs  only  a  few  times  throughout  a  flight  (and  can  be 
aperiodic  and  data-driven),  noisy  input  data  may  be  filtered  for  a  relatively 
long  time  before  use,  reducing  the  effect  of  tracking  noise.  Also,  as  the  OPTG 
technique  is  of  the  two-point  boundary-value  class,  an  intercept  goal  is  always 
before  the  projectile.  Thus,  data  drop-outs  or  loss  of  link  with  host  system 
("maximum  noise")  can  be  dealt  with  in  a  straightforward  manner:  stay  on 
course  to  the  last  predicted  intercept  point. 

7.4  Blended  OPTG /Baseline  Guidance 

To  compare  the  guidance  techniques  in  closed-loop  (say,  against 
maneuvering  targets),  it  is  necessary  to  modify  the  prototype  OPTG  guidance 
by  switching  to  its  competitor  (the  baseline  method)  in  the  end-game  (as 
entirely  OPTG  closed-loop  guidance  synthesis  is  beyond  the  scope  of  this 
Phase  I  feasibility  study).  For  this  reason,  the  closed-loop  comparison  results, 
especially  in  squib  conservation,  are  slightly  less  dramatic  than  would  be 
possible  with  a  pure  OPTG  method  (such  as  that  proposed  for  Phase  11)  which 
employs  target  feedback  information;  nevertheless,  the  improvements  are 
significant. 


In  this  "blended"  guidance  law,  the  OPTG  system  of  Figure  3.b  becomes 
the  initializing  component,  and  the  baseline  law  (of  Figure  2)  is  used  in  the 
end-game.  Switching  to  the  baseline  guidance  method  in  the  final  second  of 
flight  provides  a  workable,  if  sub-optimal,  way  to  follow  the  target  for 
preliminary  evaluation  purposes.  The  performance  of  a  complete  OPTG 
system  can  clearly  be  lower-bounded  by  this  blended  version;  still,  as  shown 
below,  even  it  outperforms  what  is  a  strong  conventional  technique:  the 
pure  baseline  guidance  method. 

The  blended  method  takes  advantage  of  the  early  OPTG  guidance 
decisions  (when  they  can  have  the  greatest  effect)  and  the  closed-loop  accuracy 
of  the  baseline  method  (which  is  good  when  the  required  divert  is  small). 
When  firing  at  stationary,  but  noisy  targets  switching  to  baseline  guidance 
from  OPTG  guidance  when  the  estimated  time-to-go,  Tau,  is  below  1.0, 
provides  the  accuracy  results  summarized  in  Figure  9  (and  analyzed  below). 
In  short,  the  accuracy  of  the  blended  system  is  slightly  better  than  OPTG  and 
much  better  than  any  baseline  case.  As  expected,  the  average  number  of 
squibs  fired  is  intermediate  to  that  of  each  of  the  techniques  when  employed 
alone  --  the  squib  use  distributions  of  which  are  depicted  in  Figure  10.  Use  of 
a  threshold  with  the  blended  method  was  found  to  affect  these  results  only 
very  slightly. 

7.5  OPTG  Guidance  Accuracy  Potential 

As  the  firing  of  a  squib  is  a  discrete  event,  the  projectile  handprint  only 
approximates  a  continuum;  actually,  a  given  trajectory  contains  N  decision 
points  separated  by  At  seconds,  and  therefore  a  maximum  of  2^  possible 
positions  at  any  given  range.*  "Perfect"  guidance  could  thus  be  defined  as  the 
ability  to  reach  all  such  possible  final  positions.  Significantly,  in  both  theory 
and  practice,  only  OPTG  guidance  can  approach  this  limit. 

While  the  accuracy  of  conventional  guidance  techniques  degrades  in 
rough  proportion  to  the  divert  (distance  from  the  ballistic  point)  required,  it 


*  For  the  given  X2f  and  At,  N  is  relatively  large  (greater  than  60)  so  the  continuum 
approximation  is  appropriate. 
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Figure  9.a:  Baseline  Targets  with  CPA  <  3  ft. 
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Figure  9.b:  Blended  Targets  with  CPA  <  3  ft. 
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Figure  9.c:  Baseline  Targets  with  CPA  <  6  ft. 
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Figure  9.d:  Blended  Targets  with  CPA  «  6  ft. 
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Figure  9.e:  Baseline  Targets  with  CPA  <  9  ft. 
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Figure  9.f:  Blended  Targets  with  CPA  <  9  ft. 


is  possible  for  an  OPTG-guided  projectile  to  be  accurate  throughout  the  entire 
envelope  of  positions  reachable  by  the  projectile;  i.e.,  to  approach  perfect 
guidance.  This  capability  is  foreshadowed,  or  hinted  at,  by  use  of  the  simple 
initializing  APN  --  as  can  be  seen  from  Figure  9,  which  compares  baseline  and 
blended  guidance  handprint  subsets  for  various  accuracy  thresholds.  For  both 
techniques,  more  and  more  targets  can  be  defined  as  "hit"  as  the  threshold  on 
the  miss  distance  allowed  rises  from  3  to  6  to  9  feet.  With  the  baseline 
method,  the  handprints  expand  outward  from  a  "tilted"  region  above  the 
ballistic  point,  into  a  diamond  shape  ~  demonstrating  that,  for  conventional 
methods,  guidance  accuracy  and  the  size  of  the  maneuver  envelope  are 
inextricably  linked. 

With  the  blended  (OPTG/ baseline)  method,  however,  virtually  the 
entire  maneuver  envelope  is  represented  at  each  level  of  miss  distance  (or 
closest  point  of  approach).  The  handprint  fills  out  and  expands  slightly  as  the 
miss  distance  threshold  rises,  but  reveals  the  basic  circular  shape  of  the 
envelope  of  potential  maneuvers  from  the  beginning.  As  is  common  with 
APNs,  the  best  points  (the  target  locations  for  which  the  model  is  most 
accurate)  are  those  best  represented  in  the  database;  here,  a  ring  with  a 
roughly  20  ft.  divert  radius  centered  around  the  ballistic  point.  Synthesis  of  a 
more  accurate  and  complete  initializing  APN  could  easily  fill  out  and  round 
out  this  ring,  leading  to  a  3-ft.-miss  OPTG  envelope  covering  about  6400  ft.^  — 
almost  four  times  larger  than  the  1700  ft.^  area  of  the  comparable  baseline 
handprint. 

Refinement  of  the  initializing  APN  and,  more  importantly,  synthesis 
of  re-initializing  APNs  (which  would  allow  the  method  to  be  purely  OPTG), 
are  sure  to  improve  the  handprints,  and  the  authors  are  confident  that 
optimization  of  the  initialization  process  is  achievable  and  will  lead  to  system 
performance,  in  the  absence  of  other  effects,  approaching  that  of  perfect 
guidance  for  a  reasonable  definition  of  acceptable  miss  distance. 

Figure  9  compares  the  "hit  distribution"  of  both  guidance  methods  for 
a  sample  of  101  low-acceleration-noise  engagements  (mean  and  standard 
deviation  =  -0.5  ft./sec.^),  and  Figure  11  presents  that  data  in  the  form  of  "^he 
cumulative  number  of  targets  hit  for  a  given  allowable  closest  point  of 


Closest  Point  of  Approach  (CPA)  (ft.) 


Figure  11 :  Cumulative  CPA  Comparison 
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approach  (CPA)  (i.e.,  the  integral  of  the  histogram^  or  density  distribution  vs. 
CPA).  Note  that  at  every  value  of  CPA  the  blended  guidance  outstrips  the 
baseline  technique;  only  a  third  of  the  baseline  trajectories  miss  by  less  than 
five  feet,  for  instance,  while  two  thirds  of  the  blended  trajectories  do.  The 
median  CPA  of  the  baseline  engagements  is  10  ft.,  compared  to  the  blended 
median  of  4  ft.  If  a  CPA  of  10  ft.  is  allowable,  over  90%  of  the  blended 
engagements  are  hits.  As  the  following  section  reveals,  such  improvements 
are  possible  for  maneuvering  targets  as  well. 

7.6  Maneuvering  Target  Results 

Against  targets  maneuvering  according  to  Eqs.  7:1  -  7:2,  with  some 
noise  on  the  measured  accelerations  (mean  =  -.25,  a  =  -.50  ft./sec.^^  the 
blended  and  baseline  systems  performed  according  to  the  handprints  in 
Figure  12.  (Note  that,  for  comparison  purposes,  the  bounds  of  the  handprint 
region  were  made  similar  to  those  of  prior  cases  by  offsetting  the  initial  target 
xi  and  X3  positions  by  25  and  30  feet  respectively.)  The  OPTG-based  (blended) 

form  remains  quite  circular  (and  is  closest  to  matching  the  targets,  which  end 
up  in  a  square  grid,  90  feet  on  a  side),  but  the  baseline  handprint  is  distorted 
further  from  the  diamond  shape  which  it  can  at  best  attain.  The  average 
divert  of  the  blended  method  was  23  percent  greater  than  that  of  the  baseline 
method  (37.18  compared  to  30.26  ft.),  and  the  average  CPA  was  cut  in  half 
(baseline  CPA  mean:  10.62  ft.,  blended:  5.26  ft.;  baseline  CPA  standard 
deviation:  8.01  ft. /sec.^,  blended:  4.86  ft./sec.^). 

Figure  13  depicts  the  handprint  subsets  for  which  the  CPA  was  less 
than  or  equal  to  six  feet.  The  circular  area  covered  by  the  blended  handprint 
(even  with  interior  gaps)  is  roughly  3.5  times  larger  than  the  smaller, 
rectangular  region  of  the  baseline  technique,  indicating  that  OPTG  guidance 
can  expand  the  region  of  weapon  effectiveness  for  a  required  closest  point  of 
approach.  The  CPA  distribution  for  all  of  the  experiments  is  shown  in  Figure 
14.  A  simplification  of  the  graphic  information,  using  the  statistics  above, 
would  be  to  summarize  that  the  blended  CPA  is  evenly  distributed,  at  a 
height  of  ten  percent,  over  a  ten  foot  region,  whereas  the  baseline 
distribution,  at  a  height  of  five  percent,  goes  out  to  twenty  feet. 
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Figure  14:  CPA  Distribution  for  Maneuvering  Case 
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Figure  15  presents  the  accuracy  information  in  another  fashion:  the  X2 
=  10,000  ft.  plane  is  shown  riddled  by  the  engagement  "bulletholes"  (or 
"projectileholes")  which  would  result  from  placing  each  target  at  the  origin. 
Here,  the  blended  guidance  delivers  a  tight  grouping  of  shots  fairly  uniformly 
distributed  around  the  "bullseye",  while  the  baseline  guidance  is  much  more 
spread  out  and  less  accurate.  The  distinctive  "X"  pattern  of  the  latter  results 
from  its  diamond-shaped  handprint;  targets  requiring  displacement  in  both 
the  X|  and  X3  directions  are  the  most  difficult  for  the  baseline  technique  to  hit, 

whereas  they  cause  no  particular  hardship  for  an  OPTG-based  guidance 
method. 

Clearly,  a  major  source  of  the  blended  guidance  accuracy  is  its  ability  to 
deliver  well,  over  the  full  region  of  weapon  operability,  the  divert  demanded. 
Figure  16  compares  the  ability  of  the  two  techniques  to  track  the  required 
divert;  the  blended  law,  it  can  be  seen,  approaches  much  better  the  ideal  line 
of  unity  slope.  In  fact,  given  targets  with  the  maneuver  capability  as  defined 
above,  one  can  conclude  that  the  OPTG-based  method  is  accurate  enough 
over  a  wide  enough  area  to  hit  (i.e.,  have  a  CPA  of  a  few  feet)  any  target  at 
which  it  is  correctly  (ballistically)  aimed  --  regardless  of  the  general  maneuver 
direction  of  the  target.  Alternately,  for  non-maneuvering  targets,  the 
authority  of  the  OPTG  projectile  is  such  to  compensate  for  4-5  mils  of  total 
aiming  error  in  azimuth  or  elevation. 

To  summarize  the  compcirison  results,  it  is  helpful  to  define  a  "ratio  of 
improvement"  as  the  baseline  CPA  of  a  given  trajectory  divided  by  the  CPA 
resulting  from  blended  guidance.  The  median  such  ratio  for  the 
maneuvering  target  case  is  1.793.  Figure  17  shows  the  statistic  collected  into 
divert  bins  of  five  foot  width,  averaged,  smoothed,  and  plotted  onto  a  log 
scale  against  required  divert.  Note  that  improvement  is  as  much  as  an  order 
of  magnitude  in  the  intermediate  region:  30  to  50  feet.  Baseline  guidance 
appears  competitive  for  small  diverts  from  ballistic  (up  to  20  feet),  and  both 
methods  begin  to  run  into  the  physical  maneuver  limits  of  the  squib-based 
weapon  system  beyond  about  55  feet.  Overall  however,  the  OPTG-based 
guidance  technique  is  shown  to  be  more  accurate  than  the  conventional 
method  at  every  region  of  divert. 


X1  (East)  CPA  (ft.) 

Figure  15.b:  Blended  CPA  Distribution 


7.7  Required  Throughput 

Each  decision  timestep.  At,  the  OPTG  guidance  system  must  update  the 
adjoint  differential  equations  described  in  Section  4.5.  Less  often  (perhaps  one 
tenth  to  one  one-hundredth  of  the  time)  an  APN  is  employed  to  reinitialize 
the  equations  due  to  new  information  about  the  target  activity.  Summing  up 
the  processor  requirements  of  these  two  components,  therefore,  provides  a 
good  estimate  of  the  throughput  required  in  an  OPTG  guidance  law 
implementation,  exclusive  of  the  tasks  of  sensor  data  filtering,  target  tracking, 
advanced  time-to-go  or  distance-to-go  calculations,  or  explicit  computation  of 
estimated  miss  distance  (if  used). 

An  efficient  form  of  the  adjoint  equations  has  been  found  to  consist  of 
156  multiplies,  91  adds,  and  7  look-up  operations  (the  latter  for  sine,  cosine, 
and  square  root  functions).  At  a  At  of  0.05  sec.,  these  254  calculations  translate 
to  5,080  operations  per  second.  The  APN  of  Figure  3.c  contains  14  multiplies 
and  6  adds  (in  network  form),  but  a  more  comprehensive  APN  might  consist 
of  up  to  50  of  each;  i.e.,  100  operations.  Even  at  a  reinitialization  interval  of 
0.50  seconds,  such  an  APN  would  only  contribute  200  operations /sec.  to  the 
required  throughput,  for  a  total  of  5,280  floating-point  operations  per  second. 
Therefore,  a  processor  capable  of  only  5.3  kFLOPS  (a  quite  reasonable  number) 
would  suffice  for  the  OPTG  guidance  calculations. 
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APPENDIX  A 


A.  APPROXIMATE  EQUATIONS  OF  TRANSLATIONAL  MOTION 
FOR  CAT  PROJECTILE 

The  three-degree-of-freedom  translational  equations  of  motion  for  the 
CAT  projectile,  treated  as  a  mass  particle,  are,  in  a  wind-axes  coordinate 
system,  using  a  flat-Earth  approximation  and  ignoring  the  effects  of  winds: 

mV  +  mg  sin  Y  +  q  S  Cd  +  T  (cos  (p  sin  -  sin  (p  sin  P^)  =  0  A;1 
mVy  cos  y  -  q  S  Cy^  -  T  sin  <p  cos  P^  =  0  A:2 

mVy  +  mg  cos  y  -  q  S  Cl^  -  T  cos  (p  cos  =  0  A:3 

in  which: 

m  =  mass  (here  assumed  constant) 

V  =  velocity  magnitude,  measured  tangential  to  instantaneous 
flight  path 

g  =  acceleration  of  gravity  (here  assumed  constant) 

y  =  flight  path  elevation  angle,  measured  positive  up  from 
horizontal  plane 

1  , 

q  =  dynamic  pressure  =  j  P 

p  =  atmosphere  density  (here  assumed  constant) 

S  =  reference  area 

Cp  =  coefficient  for  component  of  aerodynamic  drag  force  acting 
along  negative  of  velocity  vector 
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coefficient  for  component  of  aerodynamic  normal  force 
acting  perpendicular  to  the  velocity  vector  in  the  vertical 

plane  containing  measured  positive  up 
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Cyw  =  coefficient  for  component  of  aerodynamic  normal  force 
acting  in  a  horizontal  direction  that  increases  \{/ 

T  =  thrust  perpendicular  to  projectile  body  it  axis  (spin  axis), 
acting  at  roll  attitude  angle  (f> 

(p  =  projectile  roll  attitude  angle,  measured  clockwise  from  the 

vertical  plane  containing  it  (as  viewed  from  rear  of 
projectile) 

=  aerodynamic  angle  of  attack  in  wind-axes  system 
Pyy  =  aerodynamic  angle  of  sideslip  in  wind-axes  system 

V  =  flight  path  heading  angle,  measured  (verticle)  axis 
( )yf  =  denotes  the  wind-axes  system 


For  small  and  pyy,  the  lift  and  side-force  coefficients  may  be 
expressed  as: 

P 

“  ^N(x  ®w  ~  V  Pw 

P 

^Yw  ”  “  ^Na  Pw  "  V" 
where: 

^Na  =  aerodynamic  normal  force  coefficient  (here  assumed 
constant) 

Cj  =  Magnus  force  coefficient  =  d  ref  (here  assumed 

constant) 


dref  =  reference  length 

P  =  rolling  rate,  measured  positive  in  RH  direction  along  projectile 

body  1?  axis,  i.e.,  positive  in  direction  of  increasing  9  (P  is  here 
assumed  constant) 

The  drag  coefficient  may  be  written  as  a  quadratic  drag  polar: 

=  Cdq  +  K  (Cl^  +  Cy^)  A;6 

where: 

Cdq  ~  coefficient  when  lift  and  side  forces  are  zero  (here  C^q  is 

assumed  constant) 

K  =  induced-drag  factor  (here  assumed  constant) 

Substituting  Eqs.  A:4  -  A:6  into  A:1  -  A:3,  one  obtains,  for  the  case  in 
which  ot^  and  are  small: 

mV  +  mg  sin  Y  +  qS  [Cq^  +  +  p^)] 

+  T  cos  9  --  P^  sin  9)  =  0  A:7 


f  P  ^ 

mV V cosy  +  qSj^CN^jjPw  +  Tsin9  =  0 

mV  Y  +  mg  cos  y  -  qS  -  C,  ^P^J  -  T  cos  9  = 


0  A:9 


Also,  for  small  and  P^: 


<Xw  =  atot«>s9 


V  =  -  «tot  sin  9 


where  is  the  total  aerodynamic  angle. 
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In  Ref.  2  it  is  reported  that  the  firing  of  a  squib  in  the  CAT  projectile 
produces  two  temporary  effects:  a  brief  reaction  force,  T,  of  duration  6t  (a  few 
milliseconds)  and  a  transient  aerodynamic  force.  When  the  transient  forces 
have  decayed,  "...the  resulting  steady-state  solution  to  the  change  in  heading 
from  lift  force  effects  is  given  by 

in  which  is  the  maximum  value  of  the  damped  oscillation,  (Oj  is  the 

XXiaX 

nutation  frequency  (the  frequency  of  the  projectile  short-period  oscillation), 
and  CO2  is  the  precession  frequency. 


But  mVA6  is  the  change  in  momentum  of  the  particle  and  is  equal  to 
the  effective  lift  force  multiplied  by  a  time  interval  At  established  by  the  time 
required  for  the  transient  aerodynamic  force  to  subside.  However,  the 
effective  lift  force  may  be  expressed  as  q  S  where  is  the  effective 

angle  of  attack  over  the  interval  At  due  to  the  thruster  impulse.  Therefore 


qSCN^aeffAt  =  mVA0 


A:13 


and 


“eff 


mVA9 

qSCN^At 


A:14 


From  Ref.  2 
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Igliexp 

'-<1/2  ^ 

V'  - 

wherein: 


A:15 


Ig  =  squib  total  impulse  (approximately  0.75  Ib.-sec.)  =  T  5t 


Ij  =  thrust  moment  arm,  measured  positive  ahead  of  C.G. 


lyy  =  pitching  moment  of  inertia  of  projectle 
=  nutational  damping  coefficient 


Combining  Eqs.  A:12,  A:14,  and  A:15: 

tteff  = 

Let 

T*  =  scaled  thrust  =  T  5t/At  «  T 

and 


«eff  = 

where 

0)  - 


A;16 


A:17 


A:18 


A:19 


For  present  purposes,  Cq  is  assumed  to  be  constant.  But  may  be  viewed  as 
the  total  aerodynamic  angle  during  the  interval  (t,  t  +  At),  and  thus  Eqs.  A:10, 
A;ll  become,  for  the  case  in  which  a  squib  is  fired  at  the  beginning  of  the 
interval: 
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Having  accounted  for  the  aerodynamic  effects  of  a  squib  firing,  and 
noting  that  the  reaction  effects  are  those  of  a  scaled  thrust  T*  acting  over  the 
entire  interval  (t,  t  +  At),  one  has  (from  Eqs.  A;7  -  A:9)  the  following  two  sets 
of  equations  of  motion. 


When  a  squib  is  not  fired 
mV  +  mg  sin  Y  +  qSC^Q  =  0 

mV  V  cos  Y  =  0 
mV  Y  +  mg  cos  Y  =  0 
When  a  souib  is  fired 


mV  +  mg  sin  Y  +  q  S  Cj 


+  Co(T*)2  [qSK  +  Cf^  Co  +  Hj]  =  0 


mV  V  cos  Y  +  T*  (q  S  C0H3  -  sin  9)  =  0 
mV  Y  +  mg  cos  y  +  T*  (q  S  C0H5  -  cos  9)  =  0 


wherein: 

H, 


=  -  Cnoi  sin  9  +  Cl  COS9 


^5  “  -  CNotCos9  -  Cj  y  Sin9 


A:22 

A:23 

A:24 


=  0 

A:  25 

A:26 
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For  present  purposes,  we  shall  employ  the  following  parameters, 
treating  them  as  constants.  In  future  work,  these  may  be  taken  as  functions  of 
the  indicated  variables: 

^  ®  ^  ®  ^YPa  ^ref  (jv) 

C3  =  k(Cn„2  +  C22)Co2 

« 

Next,  let  N  be  a  variable  that  signifies  if  a  thruster  is  fired  within  (i.e.,  at 

*  • 

the  beginning  of)  a  solution  interval  (t,  t  +  At).  Then  N  is  binary,  N  =  0  in- 
dicating  no  firing,  N  =  1  indicating  that  a  firing  occurs.  Now,  the  equations 
of  motion  (A:22  -  A:27)  become: 

mV  +  mg  sin  Y  +  qSC^Q  +  (T*)^(qSC3  CoHj)N  =  0  A:33 

mV  V  cos  7  +  T*  (q  S  C0H3  -  sin  (p)  N  =  0  A:34 

mV  y  +  mg  cos  y  +  T*  (q  S  C0H5  -  cos  (p)  N  =  0  A:35 


A:31 

A:32 


In  Eq.  A:33,  N  is  not  squared  because  its  influence  on  the  thrust  terms  is  the 
same  whether  used  as  N  or  N^. 


APPENDIX  B 


B.  COEFnCIENTS  IN  PROJECTILE  DRAG  POLAR 

The  drag  polar  of  a  projectile  is  the  series  expansion  of  its  wind-axes 
drag  coefficient  (C15)  as  a  function  of  its  wind-axes  lift  coefficient  (Cl): 

^D  “  ^Do  I^lI  ^d/^L  ^L^^d/^L^  ■*■  ••• 

Typically: 

Cd  =  ^Do  ^Da2 

Cl  =  Ci^a  +  Clo^3  a3  + ...  B:3 

where  a  is  here  used  to  represent  the  effective  total  aerodynamic  angle  of 
attack. 

Wind-tuimel  data  are  customarily  expressed  in  the  form  of  projectile 
body-axes  force  coefficients.  As  shown  by  Hutchings*  partial  derivatives  in 
the  wind-axes  system  may  be  written  in  terms  of  partial  derivatives  in  body 
axes,  viz.: 

Cdo  =  ^Xo 

CDa2  =  ^Na  '  I^Xq  +  ^Xa2 
^La  =  ^Na  '  ^Xq 

CLa3  =  ^Na3  "  ^Xa2  "  ^^Na 

where  the  subscripts  X  and  N  denote  axial  (rearward)  and  normal  force 
components,  respectively. 


Hutchings,  Thomas  D.,  notes  dated  September  13, 1983 


Now: 


aco/acL  = 


acp/acL^  = 


^p/^g  _  ^Da_  Q 
3Ci5/3a^  ^Dg2 

dCi}/da^  "  aq^/^ 


But,  from  Eq.  B:3: 

Cl^  =  (Cl<»)2«2  +  - 

whence 


3Cl2/3o2  =  (Ci^)2 
and,  therefore 


B;8 


B:9 


B:10 


B:ll 


0  ^Dg2 

acD/acL^  =  ^ 

Defining 


K  s  dC^/dCi}  = 


induced-drag  factor 


and  introducing  the  quadratic-form  drag  polar 

Cp  =  ^Do 

one  then  has  C^q  determined  by  Eq.  B:4  and 

^Ng  -  I^xo  +  Cx„2 

K  = - — T - 


B:12 


B:13 


B:14 


B:15 


APPENDIX  C 

C.  PROJECTILE  AERODYMAMIC  COEPHCIENT  VALUES 

The  following  aerodynamic  coefficient  values  have  been  calculated 
from  ARDEC  data  provided  December  10,  1987: 


M 

o 

Q 

U 

K 

^Na' 

0.5 

0.467 

1.335 

5.03 

0.9 

0.514 

1.362 

5.03 

1.1 

0.582 

1.569 

5.03 

1.2 

0.631 

1.768 

5.03 

1.5 

0.572 

1.562 

5.03 

2.0 

0.500 

1.110 

5.55 

2.5 

0.444 

0.940 

5.61 

3.0 

0.394 

0.942 

5.55 

3.5 

0.336 

0.858 

5.35 

4.0 

0.308 

0.812 

5.01 

The  above  coefficient  values  are  based  upon  a  reference  area  (S)  of  0.12151  ft.^. 
The  induced-drag  factor  (K)  decreases  with  a,  indicating  that  the  drag  polar  for 
this  projectile  is  not  a  pure  parabola.  Because  (Eq.  A:  18)  is  a  small  fraction 
of  1  deg.,  we  have  used  ARDEC  data  for  a  =  0  and  a  =  1  deg.,  rather  than 
ARDEC  data  for  a  =  0  and  a  =  2  deg.,  in  computing  K. 

Figure  C.l  shows  the  dependencies  of  K,  Cq^,  and  Cjsj^  on  Mach 
number. 
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APPENDIX  D 


D.  ALTERNATIVE  ESTIMATE  OF  TIME  TO  GO* 


It  has  been  noted  that 


PSCd.(M)V^/2 

P  =  4m 


is  approximately  constant  over  the  Mach  number  range  of  interest  (1.2  <  M 
<  3.5),  and  it  can  therefore  be  shown  that 


X  =  - ' 

V  -  P 

where  Rg^  is  the  distance  to  go  to  reach  the  target. 

I 

For  the  projectile  described  in  Appendix  C,  with  m  =  0.966  slug,  p 
calculated  using  sea  level  standard  atmosphere  density  is  1.695  x  10'^.  This 
gives  close  agreement  with  the  actual  x  values  during  OPTG  guidance. 


Suggested  by  Albert  Rahe  of  ARDEC,  December  10,  1987. 
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E  SOURCE  LISTING  OF  SIMULATION  PROGRAM 


The  critical  components  of  the  3DOF  simulation  are  listed  here.  As  the 
program  is  in  the  "C"  language,  the  file  in  which  a  routine  appears  is 
significant,  so  a  brief  directory  of  routines  and  their  locations  is  provided. 
The  "makefile",  appearing  last  (page  115),  reveals  the  file  dependencies  and 
automatically  governs  program  re-compilation. 


Eilfi 

Routine 

Page 

main.c 

main 

89 

getdir 

91 

publish 

91 

record 

92 

fly.c 

conds  allow 

93 

fly 

94 

finalize 

96 

deriv 

97 

set_st 

99 

pull_st 

99 

tconst 

99 

advise.c 

more_trials 

100 

advise 

102 

cascade 

103 

dyn.c 

dyn 

104 

quad.c 

quad 

107 

rng.c 

gauss 

108 

random 

108 

InclHdeJile 

Page 

standard.h 

109 

integrate.h 

110 

system.h 

111 

init.h 

112 

fandb.h 

113 

bounds.h 

114 

s 


IfiV 


% 


R 


I 

ff.* 


i include  "standard. h"  /* 

/*  (stdio.h,  math  £ns,  deg,  rad,  nl)  */ 
•include  "integrate. h"  /* 

/*  {STATE_VAR,  N_STATES,  StO,  LOOP_St(), 
•include  <interfacel . h>  /* 

•include  <Getf i lePkg . h>  /* 

•include  <screen.h>  /* 


standard  functions,  consts 


integration  structure 
FORK  )  ) 

interface  strct  DESCRIPTION 
GETFILE,  FP,  GFILE,  gclose 
wide,  unwide  fns 


/*  Global  variables  used  in  dynamics 

/*  state  variables 

double  phi; 

double  x[4] ; 

double  V,  gamma,  psi; 

double  lam [13]; 

double  xT[4],  vT[41; 


rotation  angle  (rad) 
position  vector  (ft)  (1-3) 
velocity  vector  (ft/s,  rad) 
lambdas  (1-12;  4-6  const) 
target  position,  velocity 


/*  pseudo-state  variables 

double  P; 

double  W[ 9 1 ; 

double  SW[4i ; 

int  N; 

double  N_dot; 
double  uT[ 4 ] ; 
double  Umean,  Usig; 


roll  rate  (deg/sec) 
guidance  par2imeters  (0-8) 
guidance  parauns,  fn{W6,7,8) 
number  of  squibs  fired 
control  variable 
target  acceleration  (E,N,H) 
T  acceleration  noise  params 


double  Range,  Tau; 
int  hunting,  S_last; 
int  squib; 
int  Abort; 


/*  range,  time  to  target  (dyn)  •/ 
/*  dt  iteration  vars  (fly)  */ 
/•  flags  squib  firing  (fly)  */ 
/*  halt  flag  */ 


iterface  vars 


/*  (fwdinfo,  backinfo;  included  after 


Larations  &  interface  include)*/ 


double  duration; 

/* 

max  flight  duration  (sec) 

*/ 

double  dtO; 

/* 

original  time  step  (sec) 

•/ 

int  debug; 

double 

ddebug ; 

/* 

debug  flag  (used  everywhere) 

*/ 

int  bounce; 

double 

dbounce; 

/* 

"bounce"  flag 

*/ 

int  G_mode; 

double 

dG_mode; 

/* 

guidance  mode  (0,1,2) 

*/ 

int  R_mode; 

double 

dR_mode; 

/* 

recording  mode  (0,1,2) 

*/ 

int  S_mode=FALSE 

;  double 

dS  mode; 

/* 

search  mode  ( 0 , 1 , 2 , 3 , 4 , 5  ) 

*/ 

/*  (S_mode  used 

before  interface  call) 

*/ 

int  wide; 

double 

dwide; 

/* 

wide  screen  print  option 

*/ 

int  p_step; 

double 

print_int; 

/* 

print  interval  (iter,  sec) 

*/ 

int  hist; 

double 

dhist ; 

/* 

store  var  histories  flag 

*/ 

V. 

double 

dN_f  ; 

/* 

final  N  (backward  only) 

*/ 

int  integ; 

double 

dinteg ; 

/* 

integration  (0;rec;  liHeun) 

*/ 

double  Nd_thresh 
•include  "init.h" 

r 

/* 

/* 

firing  threshold 
interface  information 

*/ 

*/ 

•"j 

•  n. 

ft 


S' 

•y  .y 
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/*  (primarily)  local  shared  variables 
FILE  *fp,  *vfp; 

GFILE  *hp; 
double  Time,  dt; 
int  iter,  trials; 
int  fwd; 

double  random(),  r_time: 

double  score; 

char  strtSCR  WIDTH]; 


*/ 

/*  record,  publish  file  ptrs  •/ 
/•  gfile  ptr:  variable  history  */ 
/*  current  time  (cap!),  step  •/ 
/*  t  iterations,  trials  */ 
/*  flight  direction  flag  •/ 
/*  random  fn,  record  time  */ 
/*  trial  performance  (grope)  */ 
/*  reply  string  */ 


main ( ) 

{ 

STATE_VAR  State IN_STATES ] ; 
int  getdir();  ~ 
int  more_tr ials ( ) ,  conds_allow( ) ; 
•define  max VARS  44 
double  bstore[maxVARS ] ; 
double  fstoreimaxVARS ] ; 
int  btrials>0,  ftrials>0; 

FILE  *ffp,  *bfp; 
double  *8torage,  *sptr: 
struct  DESCRIPTION  *info,  *dptr; 
int  wscreen,  tmp; 

•  define  BANNER(h)  printf  ("\n«« 
•define  SPR(h)  8printf(str,  ((h) 
•define  ENDnote(h)  if  (lAbort)  ( 
/*  . 


/* 

129-00  simulation 

program 

*/ 

/* 

state  variables  (integrated) 

*/ 

/* 

flight  direction 

fn 

*/ 

/* 

control  routines 

V 

/* 

max  •  interface  variables 

*/ 

/* 

bJcwd  interface  var  storage 

*/ 

/* 

fwd  interface  var  storage 

*/ 

/* 

•  trials  flown  in 

each  dir 

*/ 

/* 

fwd,  bkwd  file  ptrs  (see  fp) 

*/ 

/* 

storage  pointers 

*/ 

/* 

ptrs  to  interface 

structure 

V 

/* 

screen  flag,  woric 

integer 

*/ 

1  %S  33ass:\n‘',  h) 

?  "Launch  %d"  :  "Impact  %d"),  trials) 
SPR(h);  BANNER(str);  record ( State ) ;  ) 
.  */ 


/*  un-wide  the  screen  (to  put  all  runs  on  common  footing)  */ 
scrUnwide;  wscreen  =  FALSE;  /*  (&  trac)r  screen  condition)  */ 

while(  getdir()  )  {  /*  get  direction  flag  (fwd)  for  series  */ 

/*  set  direction-dependent  interface  variables  */ 
if  (fwd)  (  storage»fstore;  info=  fwdinfo;  trials>:ftrials;  ) 

else  {  storage=bstore;  inf o=baclcinf o;  tr ials=btr ials ;  ) 

printf(  "(•trials  flown  in  that  direction  =  %d)\n",  trials  ); 

/*  open  data  base  file  only  if  going  in  that  direction  •/ 
if  ((trials)  if  (fwd)  ffp  =  fopen(  "impact. db",  "w"  ); 

else  bfp  =  fopen(  "launch, db",  "w"  ); 

fp  >■  (fwd  ?  ffp  !  bfp);  /*  set  file  pointer  */ 

while(  more_trials(  fwd,  trials,  info  )  )  { 

/*  store  interface  variables  for  recall  */ 

for  (  sptr>storage,  dptr^info;  dptr->var ! =NULL; 

*(sptr++)  =  *(  (dptr++)->var  )  )  ; 

/*  translate  interface  variables  into  internal  units  */ 
tmp  =  ( int )  (pr  int_int/dtO  ) :  p_step  =  maxd,  tmp  )  ; 

G_mode  “  ( int )dG_mode;  wide  =  (int)dwide;  debug  =  (int)ddebug; 
R_mode  =  (int)dR_mode:  hist  *  (int)dhist;  bounce  =  {int)dbounce 
S_mode  *■  (int)dS_mode;  integ=  (int)dinteg; 


/*  translate  or  initialize/f inaltze  */ 

if  (fwd)  {  psi  DEGtoRAD;  phi  *=  DEGtoRAD;  N  =  0;  ) 

else  {  psi  =  0.0;  phi  =  0.0;  N  =  (int)dN_f;  ) 

/*  Change  screen  width  if  wrong;  open  history  file  if  asked  */ 
if  (  wide  £6  Iwscreen)  (  scrWide;  wscreen<=TRUE;  ) 

else  if  ( Iwide  66  wscreen)  {  scrUnwide;  wscreen=FALSE;  ) 
if  (hist)  { 

if  (!(vfp  *  GETFILEC  hp,  "flight  data  file:",  "var.hist",  "w"))  ) 
(  pr intf ( "?Couldn' t  open  var  history  file.\n");  hist=FALSE:  ) 

) 


/*  Reset 

time 

variables 

1,  position. 

counters,  flags,  constants 

*/ 

if  (fwd) 

{  dt 

-  dtO; 

Time  «  0.0; 

FOR3(  x[k]  -  0.0;  ) 

) 

else 

{  dt 

»  -dtO; 

Time  =  duration; 

FOR3(  x[k]  «  xT[k) ;  ) 

) 

tconst ( ) ; 

/* 

squib 

constants  (fn(P,dt0)) 

*/ 

iter  * 

0; 

Abort  » 

FALSE; 

/* 

start  of  iterations 

*/ 

S  last  •• 

1; 

hunting  « 

FALSE; 

/* 

don't  hunt  'til  end 

*/ 

N_dot  -  0. 

0; 

squib  > 

FALSE; 

/* 

default:  no  firing 

*/ 

if  (R_mode 

«=1) 

r_time  = 

■■  random  (  2* 

dt,1.0 

);  /*  recording  1 

*/ 

i f  ( fwd )  ( 

/* 

set  direction  weights 

*/ 

SW[11  =  -W[61*sin(psi)  -  W(7]*cos(psi ) ;  /•  for  lamdotlO*/ 

SW(2i  =  -W(6)*cos(psi)  +  W[71*sin(psi ) ;  /*  for  laundotll*/ 

SW(31  =  W(8);  /*  for  lamdotl2*/ 

} 

else  FOR3(  SW[k]  -  0.0;  )  /*xx(not  ready  to  use  SW  bkwd )  */ 

/*  Establish  integration  structure  states  and  derivatives  */ 

if  (!fwd)  finalizeC);  /*  finalize  V,  lambdas  */ 

/*  sequester  states;  calculate  (6  seq.)  resulting  derivatives  */ 
set_st < state ) ;  deriv(state) ;  /*  (also  N_dot, squib, Range, Tau  )  */ 

ENDnote(  fwd  ) 

while  (  conds_allow( )  )  fly(  state  ); 

ENDnote ( I f wd  T  publish(); 


/*  evaluate  trial  (for  grope)  »/ 

if  (fwd)  score  ■  Range;  /*  minimize  final  range  •/ 

else  score  =  x(2];  /•  maximize  downrange  (N<Nf)  */ 


if  (bounce>s0)  ( 

/* 

bouncing  flight! 

V 

fwd  =  !fwd;  tmp=G_mode; 

/* 

reverse  6  save  flag 

V 

G_mode  =  bounce; 

/* 

use  bounce  guidance 

*/ 

sprintf(str,  "Bounce  guidance  mode:  %d",  G_mode  );  BANNER(  str  ); 


iter 

-  0; 

Abort 

-  FALSE; 

/* 

start 

of 

iterations 

*/ 

S_last 

-  1; 

hunting 

-  FALSE; 

/* 

don '  t 

hunt  'til  end 

*/ 

N_dot  ■ 

0.0; 

squib 

=  FALSE; 

/* 

default: 

no  firing 

*/ 

/*  XX (these  3  init  lines  are  repeats;  could  become  a  func)  */ 

/*  stutter-step  to  get  back  on  dtO  multiples  */ 

dt  “  ( int )  (Time/dtO )  *dt0  -  Time;  if  (fwd)  dt  +■  dtO;  fly(  state 

dt  *  (fwd  ?  dtO  :  -dtO);  /*  now,  regular  steps  */ 

while  (  conds_allow( )  )  fly(  state  ); 

sprintf(str,  "Bounce  %d",  trials  );  BANNER(  str  ); 

record (  state  );  publish();  /*  end  of  bounced  trial  */ 

G_mode"tmp;  fwd  =  !fwd;  /*  recover  flag  6  reverse  */ 


trials  “  (fwd  ?  ++ftrials  ;  ++btrials); 
fflushifp);  if  (hist)  gcloso(hp); 


/*  trial  over  •/ 
/•  flush  files  */ 


/*  recall  original  values  for  interface  •/ 

for  (  sptr>storage>  dptr«infoj  dptr->var ! *NU1.L; 

* ( (dptr++)->var )  =  *{sptr++)  )  ; 

)  /*  (more  trials)  */ 

)  /*  (direction)  */ 

if  (ftrials)  fclose(  ffp  );  /*  close  at  run  end 

if  (btrials)  fclose(  bfp  ); 

if  (wscreen)  {  scrUnwide;  wscreen=PALSE;  )  /•  unwide  screen 
)  /*  (main)  */ 


getdir()  /*  return  flight  direction  */ 

{ 

printf(  "\nDirection  to  fly  (forward,  backward,  quit) 
if  (strt01=='f '  )  (  fwd  =  TRUE;  return(  TRUE  ) 

else  if  (8tr[0]=>*'b' )  (  fwd  =  FALSE;  returnt  TRUE  ) 

else  if  (str[0i=='q' )  {  return(  FALSE  ) 

else  {  printf(  "(Need  f,  b,  or  q  as  1st  letter)\n"  ); 


gets (  str  ); 


return(  getdir()  ); 


publish!)  /*  screen  output  routine  */ 

{ 

/*  print  column  headers  when  debugging,  screen  is  full,  or  done  */ 

if  (  debug  ||  Kiter  %  ( int ) (p_step* (SCR_HElGHT-4 ) ) )  ||  iconds  allow!)  )  { 
printf(  “\n%4s  %4s  %7s  %3s  %4s  %5s  %5s  %5s  %5s  %5s  %5s  %6s"%6s", 

"Time",  "  Tau",  "  N  dot  ",  "Phi",  "V  ",  "Gamma",  "  Psi  ", 

"xl(E)",  "x2(N)",  "x3(H)",  "Range",  "  Lam  1",  "  Lam  2");  /*  77  */ 
if  (wide)  printf("  %6s  %6s  %6s  %6s  %6s  %6s  %6s",  “  /*  49  */ 

"Lam_3",  "Lam_7",  "Lani_a",  "Lam_9",  "LamlO",  "Lamll  ",  "Laml2");/*»126*/ 

^  nl;  /*  (be  sure  to  match  column  headings  with  variables!)  */ 

/"  screen  output:  angles  in  degrees  */ 

printf(  "%4.2f  %4.2f  %7.1f%ls  %3.0f  %4.0f  %5.2f  %5.2f",  /*  39+  •/ 

Time,  Tau,  N_dot, (squib?"*": "  "),  deg(phi),  V,  deg(gamma),  deg(psi)  ); 
printf(  "  %5.0f  %5.0f  %5.0f  %5.0f  %6.0f  %6.0f",  /*  38+  */ 

x(l),  x(2],  x(3].  Range,  lamll),  laml21  );  /*=  77+  */ 

if  (wide)  printf(  "  %6.0f  »6.0f  %6.0f  %6.0f  %6.0f  «6.0f  «6.0f”, 

lam[3],  lam[7] , lam[e] ,lam[9] ,  lam( 10 ) , laml 11 1 , lam ( 12 1  ) ;  /*  49+  */ 

nl;  /*=126+  */ 

if  (hist)  {  /*  file  output:  angles  in  radians  */ 

fprintf!  vfp,  "  %.31f  %.31f  %.31f  %.31f  %.31f  %.31f  %.31f". 

Time, Tau,  N_dot,phi,  V, gamma, psi  ); 
fprintff  vfp,  "  «.31f  %.31f  %.31f  ».31f\n",  x(l ] ,xl 2 ) ,x[ 3 J ,  Range  >; 
if  (G_mode*"2)  (  /•  OPTG  guidance;  show  leunbda's  •/ 

fprintf!  vfp,  "  %.31f  %.31f  %.31f  %.31f  %.31f  %.31f", 
lamll] , lam( 2 ] , lam( 3] ,  lam(7] , lam(8] , lam[9]  ); 
fprintf!  vfp,  "  «.31f  %.31f  %.31f\n-,  lamllO] , lam[ll] , lam[12]  )  ; 

)  fprintf!  vfp,  "\n"); 

)  /*  (hist:  use  tabs,  not  spaces  if  cricket;  dbp  wants  decimal  pts )  */ 


record(  state  )  /*  file  conditions  */ 

STATE_VAR  stated;  /•  state  variables  •/ 

{ 

/*  Recording  mode:  0;  full  states  (plus)  */ 

/*  (R_mode)  1:  throughout  flight  •/ 

/*  ~  2:  simple  study  (a_i  vs.  b_f)  V 

/*  (overridden  by  S_mode  when  "firtree"  generating)  */ 

int  start  ■  liter;  /*  start  of  flight  V 

static  int  Pirster*TRUE;  /*  first  of  2  trials  */ 

tdefine  fnl  fprintf(fp,  "\n")  /*  file  nl  macro  V 

#define  PRINT  St  (  LOOP  St(  fprintf(fp,"  %g",  St(lt));  )  fnl;  ) 

/*  . . . . . . V 

if  (S_mode”5)  {  /*  output  firtree  database  */ 

if  (trials  l-O)  {  /*  (after  initial  trial)  */ 

if  (Firstar)  /*  1st  flight  */ 
if  (start)  PRINT_St 

else  fprintf(fp,"  %g  %g\n",  xtl]>x(2]  );  /*  1st  impact  */ 


else  if  (istart)  fprintf(fp,"  %g  %g\n\n",  x(l),x[21  );  /*  2nd  •/ 

if  (start)  Firster  -  iFirster;  /*  toggle  */ 

)  /*  trials  V 

)  /*  firtree  S_mode  */ 

else  if  (fwd  (R  mode»2))  {  /*  simple  study  */ 

if  (start)  fprlntf (fp,"%3d.  %8.21f  %8.21f  %8.21f  ", 

trials,  laffi[4],  lam[5],  lamtS]  ); 
else  fprintf (fp, "tad.  %8.11f  %8.11f\n", 

N,  xll],  x[31  ); 

)  /*  (decimal  points  are  for  the  sake  of  dbprep)  */ 

else  {  /*  normal  R_inode{s)  »/ 

if  (start)  fprintflfpT  "  %g  %g  %g  %g  %g  %g  %g  %g  %g\n", 

W(0],Wtl],W(2J,Wl31,W(4),Wl5),W(61,Wl7),W[8)  ); 

PRINT  St 

if  (istart)  fprintf(fp,  "  %g  %g  %d  %g  *d  %d\n\n", 

Time,  N_dot,  N,  Range,  Abort,  trials  ); 

)  /•  S_mode  */  ~ 

if  ( lconds_allow( ) )  {  /*  trial  done  */ 

printf(  "\t(Squib8  fired:  %2d;  Distance  from  origin:  %6.1f)\n", 

(fwd  ?  N  :  (int)dN_f-N),  R5SS (x [ 1 ] , x [ 2 ) , x I  3 )  )  ); 
printf(  "\t(Target  position;  %8.2f,  %8.2f,  %8.2f)\n", 

XT[1],  xT(2],  xT(31  ); 

printf(  "NtlTarg.-proj.  pos:  %8.2f,  %8.2f,  %8.2f)\n", 
xTm-x(ll,  xT(2)-x(2),  xT(3]-x[3)  ); 

) 

fflush(  fp  ) ;  /*  close  observation  •/ 

)  /*  record  */ 


I 

I 

I 


B 


I 


I 


/*  fly.c 


‘/ 


tinclude  "standard. h” 

/*  (io  S  math  fns,  EPS) 
•include  "integrate. h" 

/*  <STATE_VAR,  N_STATES, 
•include  "system. h" 


/*  (system:  S,m;  environment:  g,rho) 
/*  Global  variables 


/*  standard  functions,  consts 

*/ 

/*  integration  structure 
F0R3(),  StO) 

/*  system-specific  constants 


*/ 

*/ 

*/ 

*/ 


extern 

double  W[9] ; 

/* 

guidance  parameters 

*/ 

extern 

double  dtO,  dt.  Time, 

duration;  /* 

time  variables 

*/ 

extern 

int  Abort,  debug,  fwd 

» 

/* 

flags 

*/ 

extern 

int  N,  squib; 

/* 

•squibs  fired,  flag 

*/ 

extern 

int  G_mode; 

/* 

guidance  mode 

*/ 

extern 

int  hunting; 

/* 

dt  iteration  (main) 

V 

extern 

double  Tau; 

/* 

T_go  estimate  (dyn) 

V 

/*  states 

*/ 

extern 

double  phi; 

/* 

rotation  angle  (rad) 

V 

extern 

double  xt4]; 

/* 

position  vector  (1-3:  ENH) 

V 

extern 

double  V,  gamma,  psi; 

/* 

velocity  vector  (ft/s,  rad) 

*/ 

extern 

double  lam[13]; 

/* 

lambdas 

(1-12;  4-6  const) 

•/ 

extern 

double  xT(4],  vT(41; 

/* 

target 

position,  velocity 

*/ 

/*  state  derivatives 

*/ 

extern 

double  P; 

/*  orientation 

*/ 

extern 

double  X  dot (4] ; 

/*  position 

*/ 

extern 

double  V  dot,  gamma  dot,  psi  dot; 

/*  velocity 

*/ 

extern 

double  lamdot[13]; 

/*  (4-6  const) 

*/ 

extern 

double  uT(4]; 

/* 

(actual)  target  accel. 

*/ 

extern 

double  uTm{4]; 

/* 

(measured)  target  accel. 

•/ 

/*  File  variables:  maneuver  constants 
double  Tstar,  CO,  C2,  C3; 


cond3_allow( ) 
{ 


/*  stop  conditions  routine 


•/ 


int  shared  =  ( lAbort  &&  (fabs(dt)  >=  0.0005)  (xl3)  >-  -0.5)); 

if  (fwd)  return  (  shared  &&  (Time<duration )  ); 


else 


return  (  shared  n  (Time>0.0 


)  ): 


g 
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/*  One  iteration  of  the  system  */ 
/»  state  variables  */ 


fly(  state  ) 

STATE_VAR  state ( ) ; 

STATE_VAR  hypotN_STATES] ;  /* 

double  T_left;  /* 

double  val;  /‘ 

int  check, passed, converging;  /* 
extern  int  S_last;  /* 

extern  int  integ;  /* 

extern  int  R_mode,  p_step,iter; 
extern  double  r_tinie,  randomC); 
•define  R  MAX  1.0  /* 


Hypothetical  (fwd)  state  •/ 

bkwd  est  of  time  remaining  */ 
variable  to  end  at  rero  •/ 
dt  iteration  case  flags  */ 

dt  iteration  var  (main)  •/ 

integration  (0:rec;  l:Heun)  */ 

/*  record  mode,  step  */ 
/*  random  obs  time,  fn  */ 
max  sec  between  recordings  */ 


•define  LIST(p,g)  if  (debug)  {printf(p);  LOOP_St ( pr in tf ( "%d : %g  ",k,q(k));)  nl;) 

/* . - . —  */ 

/*  output  when  printstep  agrees  6  not  aborted  or  hunting  */ 
if  (l(  (iter  t  p  step)  ||  Abort  ||  hunting))  publishO; 
if  ((R_mode--l)  Is  (r_time<-dtO*iter ) )  /•  time  to  record  •/ 

(  record (  state  );  r  time  +=  random(2*dt0,  R_MAX);  ) 

LIST(  "statesi  ",  St  )  t!lST(  "dots:  ",  St_DOT  )~  /*  dbg  •/ 


check  -  TRUE;  /*  reduce  time_step  magnitude  near  end  */ 

if  (fwd  &&  (R_model*‘2 ) )  /*  goal;  time_to_go==0  •/ 

{  val  •  T_left  -  Tau;  )  ~ 

else  if  (fwd~6i  (R  mode*»2))  /*  goal;  x2““xT2  */ 

{  val  -  xt2);  T_left  -  (xT(2]-xt2J )/x_dot[2) ;  ) 
else  if  (Ifwd  a  (gamma  >  0.0))  /*  goal;  V»=V  muzzle  */ 

(  val  -  V;  T_left  =  (3821.5  -  V)/V_dot;  )" 
else  check  ■  PALSET  /*  (impact  half  of  bkwd  flight)  •/ 

/*  (note;  these  equations  depend  on  roughly  parabolic  paths)  •/ 

if  (check)  {  /*  look  for  dt-reducing  cases  */ 

passed  ••  (signum(val)  !■  S_last);  /•  goal  passed  •/ 

converging  »  (fabs(T_left)  <  fabs(dt));  /"step  shrinking*/ 

if  (hunting  | |  converging  | |  passed)  ( 


hunting  -  TRUE; 

/*  (dt 

reduction  guaranteed) 

•/ 

if  (converging) 

dt 

-  T  left; 

/*  trust  prediction 

*/ 

else  if  (passed) 

dt 

*=  -0.5; 

/•  go  back  half 

•/ 

else 

dt 

*-  0.5; 

/*  go  forward  half 

*/ 

last  ■  signum(val) 

t 

/*  (save  for  next  pass) 

*/ 

V 

S’ 


I 


V 

S* 

l.a 


I 


/*  Update  time  stamp  and  integrate  effective  Ndot  */ 
iter++;  Time  +-  dt;  if  (squib  (hunting)  if  (fwd)  N++;  else  N — ; 
if  (integ«“0)  {  /*  Integrate  states  using  rectangular  method  */ 

LOOP  St(  St(k)  +■»  dt*St  DOT(k);  )  deriv(state)  ; 


else  { 


/*  Integrate  states  using  Heun's  method  (1  fwd  est.) 


pjKl 

I 


:: 


/*  get  hypothetical  state  through  rectangular  integration  */ 

LOOP_St<  hypo[lc].u  *  St()c)  +  dt*St_DOT(  k  ) ;  )  der  i  v  ( hypo )  ; 

/*  use  average  of  forward  and  current  derivatives  as  slope  */ 

LOOP_St(  St(k)  +»  0.5*dt*(St_DOT(k)  +  hypofkl.up);  )  der iv (state); 


/*  rstates  fi  dots  up-to-date  in  both  forms  after  LOOP_St->der iv  pairs*/ 
)  /*  (fly)  */ 


K# 

W'**. 

?Ivlv 

m 


W 


W: 

M 


finalizeO  /*  set  final  V,  lambda  values  */ 

K 

double  CD0»  CNa,  K;  /*  Mach  "cons'. ts“*/ 

double  qS;  /*  dyn  pressure  */ 

double  atmp>btmp,ctmp,dtmp;  /*  qS  work  vars  */ 

double  Qa,Qb,Qc,  Mroot,  W;  /*  qS  work  vars  */ 

int  neg;  /*  negative  discriminant  flag  */ 

tdefine  Vmin  100.0  /*  minimum  V  f  */ 

/* - - - */ 

/*  lambda  5  and  6  are  constant  (and  from  interface)  •/ 
FOR3(  lainIk+6]  «  lam[k+91  »  0.0;  )  /*  7->12  are  0  */ 


/*  solve  for  qS_£  to  get  Vf  •/ 

/*  XXX  (temporarily  use  Mach=1.75  values  for  CNa,K,CD0)  */ 

CNa»5.29;  K-1.336;  CDO-0.536;  /*  (see  tables  in  dyn)  */ 

atmp  ■  (lam[S]*CNa  -  lam[ 6 1 *C2 ) *C0/Tstar ;  btmp  =  CDO/m; 

ctmp  <■  lam[  5  ]  *g*cos  (gamma)  -  W[4)*N*N;  dtmp  =  g*sin  (gcunma) ; 

dbg(  "tas  %g  tb;  %g  tc;  %g  td ;  %g\n“,  atmp, btmp, ctmp, dtmp  ); 

C3  »  K*(  CNa*CNa  +  C2*C2  )*C0*C0;  /*  new  firing  constant  */ 

Qa  >  atmp*btmp: 

Qb  «  atmp*dtmp  +  lamI5]*btmp  +  C3*ctmp/Tstar ; 

Qc  ■  lam(51*dtmp  +  C0*ctmp; 

dbg(  "Qa:  %g  Qb;  %g  Qc;  %g\n",  Qa,Qb,0c  ); 

quad(  Qa,Qb,Qc,  aqS,  tMroot,  ineg  ) ;  /*  (use  plus  root)  */ 

VV  ■  2.0*qS/(  rho*S  ) ;  /*  (qS  =  0.5*rho*V»V)  */ 

if  (neg  ||  (W  <  Vmin*Vmin))  { 

Abort  -  TRUE; 

printf(  "— >  qS_f  problem  (V*V  »  %g);  try  again. \n\n",  W  ); 

V  -  qS  ■  EPS;  /*  (just  to  finish  out  iter)  */ 

) 

else  V  ■  sqrt (  VV  ) ; 


/•  solve  for  lambda  4  (  -  -(lam5»Q5  +  Iam6*03)/Q0  )  */ 

lam[4]  -  (  qS*C0*(  lam(51*CNa  -  lamt6)»C2  )  +  lam(5]*Tstar  ) 

/(  qS*C3  +  Tstar*C0  ) ; 
printf(  "Lambda_4f:  %g\n",  laml4)  ); 


/*  now  assign  lambda's  1-3  */ 

lam[l]  >  -lam[4]/m; 
lam[2]  ■  -leun[6]/m; 
lam[3]  ■  -lamisj/m; 


)  /*  (finalize)  */ 


V 

*/ 


deriv(  state  )  /*  obtain  state  derivatives 

STATE_VAR  state[];  /*  state  variables 

{ 

extern  double  N_dot,  Nd_thresh;  /*  firing  variable,  threshold  */ 

extern  double  CDO,  CNa,  K;  /*  Mach-dependent  variables  */ 

double  H3,H5,  Q0,Q3,Q5;  /*  worlt  vars  */ 

double  s_gaiiuna,c_ganuna;  /*  worJc  vars  */ 

double  qS,  8_phi,c  phi;  /*  work  vars  */ 

double  CHH,  AO,  aT,  B;  /*  OPTG  vars  */ 

double  Uc[4],  Ndcjp8i,Ndc_gam;  /*  baseline  */ 

double  d_V_dot,d_gainina_dot,d_psi_dot;  /*  delta  dots  */ 

double  L26,  L35;  /*  delta  L  vars*/ 

•define  Nmax  4000  /*  maximum  f squibs  (xx  no  lim)  */ 

/* . . . */ 

pull_st(  state  );  dynO;  /•  pre-squib  state  derivatives  */ 

/*  (xx:  could  put  pull_st  in  dyn  again  along  w/  integrate. h)  */ 

s_phi  »  sin  (phi);  s_g^unma  =  s  in  (gamma); 
c_phi  »  cos(phi);  c_gamma  =  cos(gamma); 

H3  »  C0*(  -CNa*s_phi  +  C2*c_phi  ); 

H5  ■  C0*(  -CNa*c_phi  -  C2*s_phi  >;  /*  (H's  include  new  CO)  */ 

C3  -  K*(  CNa*CNa  +  C2*C2  )*C0*C0;  /*  new  firing  constant  */ 

/*  (C3  =  K*SS{H3,H5) )  */ 

qS  -  0.5*rho*S*V*V; 

QO  =  qS*C3  +  Tstar*C0; 

Q3  =  qS*H3  -  Tstar*s_phi ; 

Q5  =  qS*H5  -  Tstar*c  phi; 

dbg(  “QO:  %g  Q3:  »g~  Q5 :  %g\n",  Q0,Q3,Q5  ); 


/*  compute  N_dot  if  guiding,  not  impacting,  and  squibs  left  */ 
if  (  G  mode  66  ((hunting  ||  (Tau— 0.0))  66  ((  fwd  &&  (N<Nmax))  || 

( Ifwd  66  (N>0  ) ) )  )  ( 

CHH  »  rho*V*S*(  lam(ll*C3  +  lam[21*H3  +  lam[3]*H5  ); 

if  (G_mode»=2)  {  /*  OPTG  guidance  */ 

/*  Calculate  N_dot;  A2*Nd_dot  +  Al*N_dot  +  A0*N  =  B  */ 

/*  (Note:  coefficient  signs  opposite  those  of  notes)  */ 

/*  (A2  =  -2*W0*Tau“rr ,  but  Nd  dot  is  n.e.  zero)  */ 

A1  «  2.0*W[0];  /*  2.0*Wi0]*rr*Tau*(rr-l)  */ 

AO  -  2.0*W(4]; 

B  lamdot[l]*Q0  +  leundot ( 2 1  *Q3  +  lamdot[3)*Q5 
+  V_dot*CHH  +  P*(  lam(2)*Q5  -  lam[3]*Q3  ); 

/*  (...whether  we  use  derivatives  with  or  w/o  N_dot)  •/ 

N_dot  -  (  B  -  A0*N  )/Al;  /*  or  (B/2  -  N»W4)/W0  */ 

dbg(  “AO:  %g  Al:  %g  B:  %g\n",  A0,A1,B  ); 

) 

else  {  /*  G_mode««»l  — >  baseline  guidance  */ 

/*  acceleration  commands  zeroing  predicted  miss*/ 

FORK  Uc(k]  «  uTm[kl  +  2 . 0*  (vT(K  ]-x_dot  [  )c  ]  +  (xT[  It  ] -x[  )c  1  )/Tau  )/Tau 
dbg(  “Uc%d;  »g  ",  k,Uc(lt]  );  )  dbg("\n"): 


if  <fabs(Q3)  <  EPS)  Ndc_psi  =  0.0; 

else  Ndc_p8l  -  m*  (Uet2)*sin(p8i)  -  Uc[l)*cos(psi) )/Q3; 
if  (fabs(Q5)  <  EPS)  Ndc_gani  -  0.0; 

else  Ndc  gam  »  m* ( (Uc[2I*cos(psi)  +  Octl)*sin(psi) )*s_gamma 

(Uc[31  +  g)*c_gamma)/05; 

N_dot  ■  Ndc_gam*fabs(c_pUi)  +  Hdc_psi*fabs(s_p)ii) ; 
dbg(  "Ndcjpai;  %g  Ndc_gam:  %g\n",  Ndc  psi,  Ndc_gam  ); 

)  /•  (G  mode)  */ 


squib  <■  (N_dot  >  Nd_thresh); 
)  /*  (G  mode,  etc. )  */ 


/*  N_dot  large  enough?  */ 


if  (squib)  {  /*  physical  effects  of  squib  (treat  N_dot  as  1)*/ 

d_V_dot  =  -QO/m;  V_dot  +=  d_V_dot; 

d~gamma  dot  «  -Q5/(m*V);  gamma_dot  +=  d_gamma_dot; 

d_psi_dot  -  -Q3/ (m*V*c_gamma ) ;  psi_dot  +=  d_psi_dot; 

/*  (gamma  !*  k*pi/2)  */ 


if  (G_mode— 2)  (  /*  update  lamdots  1-3  */ 

L26  >  m*lam[21  -*■  lam[6]; 

L35  -  m*lamI31  +  lam[5];  /*  (LI 4  not  needed)  */ 

lamdotdl  +=  (  CHH  +  L35*d_gamma_dot  +  L26*d_ps i_dot*c_gamma  )/m! 


lamdot(2]  +=  L26*(  V*d_gamma_dot*s_gamma  -  d_V_dot*c_gamma  ) 

/(  m*V*c_gamma  ); 

lamdot[3]  +=  (  -L35*d  V  dot  -  L26*V*d_psi_dot*s_gam)na  >/(  in*V  )  ; 
)  /*  (G  mode  2)  */ 


)  /*  (squib) 

*/ 

/* 

sequester  the  state 

derivatives  in  the 

structure 

St 

DOT(  7)  - 

lamdot[ 

11 ; 

St_ 

DOT( 

0) 

m 

P; 

St 

■dOT(  8)  = 

leundot[ 

2) ; 

St; 

;dOT(  9)  = 

l2undot[ 

3]  ; 

St 

DOT( 

1) 

B 

x_dot [ 1 ] ; 

St" 

■dot( 

2) 

B 

x_dot [21; 

St 

DOT(IO)  = 

lamdotl 

71; 

St" 

DOT( 

3) 

m 

x_dot(3] ; 

St* 

DOT(ll)  - 

lamdotl 

8); 

St; 

_DOT(12)  = 

lamdotL 

9); 

St 

D0T( 

4) 

s 

V_dot ; 

St" 

■dot( 

5  ) 

m 

gamma~dot ; 

St 

DOT(13)  ■ 

lamdotl 10 ]  ; 

sti 

■dot( 

6) 

m 

psi~dot; 

St 

DOT(14)  “ 

leundotill  ] ; 

St_ 

.DOT(15)  - 

l^undot[12  ] ; 

St 

DOT (16) 

m 

vTtl];  St 

DOT(19) 

»  uT(l] ; 

St" 

'dOT(17) 

m 

vT[2];  St 

■dOT(20) 

-  uT[2]; 

St" 

'DOTdS) 

m 

vTt3];  St* 

’dot  (21) 

-  uT(3J ; 

)  /*  (deriv)  */ 


^^VWliU1^^^t^^U^VV^^W»Vlrpl^v1ru  wv  WVTTU  mj  ITU  m.  i»w  trjmr^wv  k'J 


set_st<  state  ) 
STATE  VAR  state[]; 


/*  set  states  in  integration  structure  for  dyn  */ 


input:  state  variables 


St(0)*  phi; 
St(l)»x[l];  St(4) 

St(2)»x[2];  St(5) 

St(3)*x[3];  St(6) 


St(4)=  V; 
St  (5  )>:geunma; 
St(6)=  psi; 


St(7)»laintll 
St(8  )>lain(2] 
St(9)»lam(3) 


St(16)-xT[l]; 
St(17)-xT[21; 
St(18)-xT[31  ; 


St(19)-vT(l) ; 
St(20)-vTl21; 
St(21)-vT[31 ; 


St(10)»lainl71 
St  (11  )=lani[  8  ) 
St(12)=lam[ 91 


St(13)»lam(101 
St(14)  =  lain[lll 
St(15 )=lam[ 121 


pull_8tt  State  ) 
STATE_VAR  8tate[l; 


pull  states  from  integration  structure 
input:  state  variables 


static  double  circle»2.0*PI; 

St(0)  -■  circle*(int) (St{0)/circle); 
if  (St(0)<0.0)  St(0)  +»  circle; 

phi  >St(0); 

x(ll»St(l);  V  -St(4);  lam(ll«St(7) 
x[21«St(2);  ganuna=St ( 5 ) ;  lamt 2 l=St (8 ) 
xI3]=St(3);  psi  -St(6);  lam(31=St(9) 


xT(ll-St(16); 

xT[2]=St(17); 

xT[31-St(18); 


vT(ll=St(19); 

vT[2]*St(20); 

vT[31»St(21); 


phi  wrap-around 
positive  version 


lam[71-St(10) 

lam[81=St(ll) 

lainI91=St(12) 


Iamtl01-St(13) 
lamtlll=St(14) 
lam[121-St(15 ) 


tconstO  /*  set  P-  and  dt-dependent  constants  */ 

{  /*  (called  in  main  prior  to  first  dyn)  */ 

Tstar  ■  0.75/dt0;  /*  I_s/dt_decision  */ 

CO  ■  6.5205e-6*Tstar/dt0;  /*  -new  squib  constant  •/ 

C2  »  P*8.1938e-4;  /*  P*0.5*CYp*d  ref/Vav  */ 


^'jr. 


/*  advise. c  */ 
tinclude  "standard. h” 

iinclude  "interface! . h"  /*  interface!  structure  DESCRIPTION  */ 

#inc!ude  <Getf i!ePltg.h>  /*  GETFILE,  FP,  GFILE  package  */ 


/*  G!obal  variables  from  main  •/ 


extern 

double  W[ 9 ] ; 

/* 

guidance  parameters 

*/ 

extern 

double  duration 

# 

/* 

mrocimum  flight  time 

*/ 

extern 

double  score; 

/* 

grope  metric 

"/ 

/*  extern  double  N_dot; 

* 

control  variable 

*/ 

/*  state  variables 

*/ 

extern 

double  phi; 

/*  rotation  angle  (rad) 

*/ 

extern 

double  x[ 4 ] ; 

/*  position  vector  (1-j;  EINH) 

*  / 

extern 

double  V,  gamma 

,  psi; 

/*  velocity  vector  (ft/s,  rad) 

*/ 

extern 

double  lam[13]; 

/*  lambdas  (1- 

12;  4-6 

*/ 

extern 

double  xT[4],  vT[4] 

f 

/*  target  position,  velocity 

■*/ 

extern 

int  S  mode; 

/* 

search  mode: 

*/ 

tdefine 

GRID  1 

/* 

explore 

regular  grid 

*/ 

•define 

GROPE  2 

/* 

use  grope  for  search 

•/ 

•define 

GAUSS  3 

/* 

explore 

gauss ian  space 

*/ 

•define 

UNIP  4 

/* 

explore 

unif.  random  space 

*/ 

•define 

FIRTREE  5 

/* 

generate  fwd-branching  db 

•/ 

•define 

loopK(h)  (int 

jj 

for  (j=0 

;  j<numK;  j++)  (h)  } 

/*  Globa!  variable  for  record  */ 

/*  double  TgoO;  *  firtree  base  time  to  go  */ 


more_trials(  fwd,  trials, 

info  )/* 

routine  governing 

trials 

*/ 

int  fwd,  trials; 

/* 

direction  flag,  •trials 

*/ 

struct  DESCRIPTION  *info; 

( 

•include  "bounds. h" 

/• 

ptr  to  interface  structure 

*/ 

/* 

Min_,  Max_,  Inc_, 

num_  (F,B) 

•/ 

double  K[max  (numF,nuinB)  ] ; 

/* 

parameter  values 

*/ 

int  numK; 

/* 

index,  •parameters 

*/ 

int  adviseO,  interface() 

;  /* 

routines 

*/ 

int  GOING; 

/* 

returned  flag 

*/ 

static  int  First  auto=TRUE;  /* 

auto,  mode  flag 

*/ 

/*  GFILE  *gp;  static  FILE 

*f  irp;  * 

firtree  file  ptr 

*/ 

/*  static  int  count-0; 

* 

branch  counter 

*/ 

/*  double  tmp; 

* 

read  place-holder 

*/ 

/*  define  DELTA  G  50 

* 

G4,G5  change 

*/ 

*/ 


if  (!S_mode)  {  /*  use  interface  */ 

char  *tit!e»"!29-00  TPBV  Guidance  Project"; 

GOING  »  interfaceC  title,  info,  (trials  ?  NORESET  :  RESET)  )  ; 
/*  ('quit'  exits)  */  /*  (reset  if  no  previous  trials)*/ 


100 


5 


is 


ra?  c'cMi  si  :  K 
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/*  XX  (firtree  mode  currently  disabled) 

else  if  (S_mode”FIRTREE)  (  *  adjust  params  through  database  * 

if  (First_auto)  {  *  open  reference  data  file  * 

First_autO  =  FALSE; 

GBTFILBC  gp,  "Base  data  filename:  ",  "base.db",  "r"); 
firp  »  FP(  gp  ) ;  *  file  pointer  * 

) 

if  (  (count%2)»0  )  (  *  first  trial  of  new  flight  * 

printf ( "Reading  conditions  for  reference  observation  %d.\n",  count/2); 

*  (reference  assumed  to  impact  at  (0,0))  * 

GOING  -(  3»«facanf (firp, "%lf  %lf  %lf",  SG4, SG5 ,STgoO )  ); 

GOING  *r«(10— fscanf  (f  irp,  "%lf  %lf  %lf  %lf  Ilf  Ilf  Ilf  Ilf  Ilf  Ilf", 
llam[6 1 ,llam[l] ,6lam[2] ,&lam[ 3] , 

IV, (gamma, 6psi,  lx[ 1] ,ix[ 2 ] ,lx[ 3 ] )  ); 


if  (GOING)  G4  +-  DELTA_G; 

else  printf (  " — >  Reading  Done  < — \n"  ); 


else  { 

GOING  -  TROE; 
G4  —  DELTA_G; 
G5  +«  DELTA  G; 


second  trial 


*  reset  G4 

*  change  G5 


printf (  “Try  G4:  Ig  and  G5 :  lg.\n",  G4,  G5  ); 
count++ ; 


(firtree)  xx  */ 


else  {  /*  adjust  parameters  automatically 

numK  ■  (fwd  ?  numF  :  numB); 


r» 

: 

i 


if  (First_auto)  {  /*  allow  user  to  verify  limits  */ 

First  auto  «  FALSE; 

print7(  "Minimums,  maximums  (and  increments)  are:\n"  );  loopKC 

if  (fwd)  printf(  "Id:  lg\tlg\t (Ig )\n",  j,  MinF[ j ) ,MaxF[ j ] , IncF[ j )  ); 
else  printf(  "Id:  lg\tlg\t (Ig )\n" ,  j,  MinB[ j] ,MaxB[ j) , IncB[ j ]  ); 

) 

/*  (Note:  vars  in  both  directions  must  match  those  in  bounds. h)  •/ 
if  (fwd)  {  /•  Target  East,  Height  */ 

K[0]-lam(4];  K(ll-lam[51;  K(2]»lam[6); 

GOING  ■  advi8e(  numK,  K,  MinF,  MaxF,  IncF  ); 
lam[4J-K[0];  lam(5]-K(ll;  lamI6J»KI2]: 

) 

else  { 

K(0]-lam(5];  K(ll-lam(6J;  K(2]-W[4); 

GOING  •  advise (  numK,  K,  MinB,  MaxB,  IncB  ); 
lam[51=K(0];  lam(6)-K(lJ;  WI41-KI2); 

) 

)  /*  (S_mode)  */ 

return (  GOING  ) ; 

)  /*  (more_tr ials )  */ 


advisee  numK,  K,  Min,  Max,  Inc  ) 

/*  advise  more  trials  about  parameter  values  */ 
int  numKj  ~  /*  #parameters  */ 

double  K(l;  /*  array  of  parameters  */ 

double  Mint],  Maxt ] :  /*  bounds  V 

double  Inctl:  /*  increments  (if  GRID)*/ 

^static  int  GOING-TRUE;  /*  FALSE  when  over  */ 

static  int  a_count-0;  /*  calls  counter  */ 

•  define  Max  Iter  1000  /*  maximuun  a_count  */ 

double  hal?:  /*  used  for  mean,  sigma*/ 

int  cascadeO,  gropeO;  double  gauss  (),  randomC); 

/* - */ 

a_count++; 

/*  distribute  samples  ...  */ 

if  (Sjnode— GRID  )  {  /*  . . .  regularly  */ 

if  (a_count“l)  loopK(  Ktjl  -  Mintjl;  ) 
else  GOING  -  cascade (  0,numK-l,  K,  Min, Max, Inc  ); 

else  if  (S_mode-=GROPE)  (  /*  . . .  with  a  purpose  */ 

GOING  -  ! grope (  score,  numK,  K,  Min,  Max  ); 

/*  (grope  counts  iterations  on  its  own)  */ 

) 

else  if  (S_mode— GAUSS)  {  /*  . . .  gauss ianly  */ 

GOING  -  (a_count  <=  Max_Iter); 

if  (GOING)  loopK(  half  -  (Maxljl  -  Mint j) )/2.0; 

do  KtjJ  -  (Minij)  +  half)  +  gau8S(  half/2.0  ); 
while  (  (Ktjl  <  Mintjl)  ||  (Ktjl  >  Max(j))  );  ) 


) 


/* 


. .  uniformly 


*/ 


else  if  (S_mode— UNIF  )  { 

GOING  -~(a_count  <-  Max_Iter); 
if  (GOING)“  loopK(  Ktjl~-  random(  Mintjl,  Maxljl  );  ) 

) 

else  ( 

printf(  " — >  Onlcnown  search  mode:  %d  (%g)  < — \n",  S  mode,  S  mode  )i 
GOING  »  FALSE; 

1 


');  loopK(  printf("%g  ",Kljl);  )  nl; 


if  (GOING)  ( 

printf(  " — >  Try  values; 

) 

else  { 

S_mode  -  FALSE;  /*  put  user  in  interface  */ 

printf(  "\n — >  (Restart  if  wishing  to  search  further)  < — \n"  ); 

) 

return (  GOING  ); 

)  /*  (advise)  •/ 


y. 


58 


I 

I I 


■l^;^^J1IVl*U»rJVUW^/>rJll^.>r»^^ulr»lf^:lr-•-.  r-jirujrv 


/*  dyn.c 


linclude  "standard. h"  /* 

/*  (stdio.h,  math  £ns,  max(),  dbg)  */ 

linclude  "system. h"  /* 

/*  (system:  S,mj  environment:  g,rho)  */ 
double  CDO,  K,  CNa;  /* 


Global  variables  for  others:  state  derivatives 


/* 

standard  functions,  consts 

*/ 

dbg ) 

*/ 

*/ 

/* 

system-specific  constants 

g,rho) 

*/ 

/* 

Mach-dependent  system  vars 

*/ 

/* 

calculated  from  tables  here 

*/ 

/"  (note:  orientation  (phi)  t  P  (phi 

dot) 

not  needed)  */ 

double 

x_dot [ 4 ] ; 

"  /* 

position  */ 

double 

V~dot,  gamma_dot,  p8i_dot; 

/* 

velocity  */ 

double 

lamdot[13] ; 

/* 

(4-6  const)  */ 

double 

uTm[ 4 ] ; 

/* 

measured  uT  "/ 

/*  Global  state  variables  from  main 

*/ 

extern 

double  x( 41 ; 

/* 

position  vector  (ft)  (1-3) 

*/ 

extern 

double  V,  gamma,  psi; 

/* 

velocity  vector  (ft/s,  rad) 

*/ 

extern 

double  lamfl3]; 

/* 

l^unbdas  (1-12;  4-6  const) 

*/ 

extern 

double  xT[41,  vT[4]; 

/* 

target  position,  velocity 

*/ 

extern  double  uT[4]; 
extern  double  Umean,  Usig; 
extern  double  gauss ( ) ; 


/*  constant  target  acc  */ 
/*  noise  mean,  sigma  */ 
/*  rng  for  acc  noise  */ 


dyn  ( ) 


given  states,  calculate  state  derivatives 


Generalized  3-DOF  Dynamics 
Copyright  1987  Barron  Assoc.,  Inc. 
John  F.  Elder  IV 


/*  Global  variables  from  main 

*/ 

extern 

double  W[9] ; 

/* 

guidance 

parameters 

(0-5) 

*/ 

extern 

double  SW[4i; 

/* 

guidance 

parameters 

(6-8) 

*/ 

extern 

double  Tau,  Range; 

/* 

Est  time, 

dist  to  target 

*/ 

extern 

int  debug,  G_mode; 

/* 

debugging 

,  guidance 

flags 

*/ 

/"  internal  variables  and  constants  */ 
double  vl4),  Uc(4];  /*  projectile  vel,  command  acc."/ 
double  d[4i,  e[41;  /•  position,  vel  differences  */ 
double  det[4];  /*  Taylor  series  diff.  w/o  U's  */ 
double  8_gamma , 8_ps i ,  I>sc,Ucsc;  /*  work  vars  */ 
double  c_gamma,c_psi,  Lcs,Uccs;  /*  work  vars  */ 
double  qS,  TT,  L26,L35,  tmp,tmp2;  /*  work  vars  •/ 


double  Mach;  /* 
int  ind;  double  prop;  /• 
static  double  Mbound[ll]B{  0.0, 


Mach  #  wrt  V_sound  at  alt  0  */ 
interpolation  variables  */ 


static  double  MbounddH^i  0.0,  0.5,  C.9,  1.1,  1.2,  1.5, 

2.0,  2.5,  3.0,  3.5,  4.0  }; 

static  double  CD0val[ll]»{  0.450,  0.467,  0.514,  0.582,  0.631,  0.572, 

0.500,  0.444,  0.394,  0.336,  0.308  ); 


$ 

I 
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static  double  Kval[lll-(  1.300,  1.335,  1.362,  1.569,  1.768,  1.562, 

1  .110,  0.940,  0.942,  0.858,  0.812  )  ; 

Static  double  CMaval[ll]-{  5.03,  5.03,  5.03,  5.03,  5.03,  5.03, 

5.55,  5.61,  5.55,  5.35,  5.01  ); 

/*  (values  at  Mach«0  added  by  hand;  could  instead  use  Glauert  approx)  */ 

/*  Glauert  (below  Mach  1)  or  Prandtl  (above  Mach  1)  approximations: 

/*  if  (Mach<Mbound [  0])  CD0_M  *  CdragO/sqrtd .0  -  Mach*Mach); 

/*  else  if  (Mach>Mbound[10] )  CD0_M  =  CdragO  +  0 . 1 415/sqr t (Mach*Mach  -  1.0); 

int  k;  /*  macro  loop  index  . */ 

•  define  FOR3(h)  for  (lc-1;  )c<«3;  k-i-t)  (h) 

/*  - V 

/*  adjust  drag  coefficient  according  to  Mach  number  profile  */ 

Mach  ■  V/1116.89;  /*  (using  sound  velocity  •  0  altitude)  */ 

/*  (max  Mach  of  4.0  ->  max  V-4457.56)  */ 

/*  find  lower  bound  index  and  proportion  of  overflow  */ 

for(  ind-0;  Mach>Mbound[  ind-t-1  ] ;  ind++); 

prop  >  (  Mach-Mbound[  indl  )/(  Mboundt  ind-fl  ]  -  Mboundtind]  ); 

CDO  ■  CDOvaKind]  +  prop*  (  CDOvaK  ind+1  ]  -  CDOvaltind)  ); 

CNa  ■  CNavallindj  +  prop* (  CNavalt ind+1 1  -  CNavaliind]  ); 

K  ■  Kvaliind]  +  prop* (  Kval(ind+ll  -  Kvaliind)  ); 

dbg(  "Mach:  %g;  CDO;  %g;  CNa:  %g;  K:  %g\n",  Mach,  CDO,  CNa,  K  ); 


qS  -  0.5*rho*V*V*S; 
s_gafflma  «  sin (gamma); 
c_gamma  ■  cos (gamma); 

/*  Geographic  rates  ' 


s_psi 

c_p8i 


sin(psi) ; 
cos (psi ) ; 


■  x_dottl] 

■  V*c_gamma*s_psi; 

/* 

"East" 

*/ 

■  x~dott2] 

■  V*c_gamma*c3psi; 

/* 

"North" 

*/ 

-  X  dot(3) 

■  V*s_gaumna; 

/* 

up 

*/ 

/*  Acceleration  vector  (sans  squib  effects)  */ 

V_dot  ■  -g*s_gamma  -  qS*CD0/m;  /*  ft/ss*/ 

gamma_dot  ■  -g*c_gamma/V;  /*  rad/s*/ 

psi_dot  -  0.0;  /*  rad/s*/ 

/*  simple  Tau  estimation  =  Range/(V  along  Range)  */ 
tmp  ■  tmp2  *  0.0;  /*  use  pos  and  vel  differences  */ 

FOR3(  d[Js]  »  xT[)c]-x()t] ;  tmp  +=  dl)c]*dU]; 

e[)s]  ■  vT[Js]-v()c] ;  tmp2  -=  d(Jc]*e[)c);  ) 

Range  -  sqrt(  tmp  ); 

if  (Range»0.0)  Tau  =  0.0;  /*  impact  */ 

else  Tau  =  (Range/(  V  -  0.001695*sqrt(V)*Range  ) )*signum(tmp2) ; 

/*  (const  ideal  for  lOK  ft,  3.605  sec  flight  with  32  squibs) 

/*  (signum  is  for  iteration  purposes)  */ 

/*  measure  (with  noise)  target  acceleration  */ 

FOR3(  uTm[)c]  ■  uTtic]  +  Umean  +  gauss  (  Usig  );  ) 


m 
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/*  Adjoint  guidance  equations 
if  (G  mode  l>  2) 


/*  baseline  or  no  guidance 


F0R3(  lamdott)c]*lamdotI>c+6]«lamdot(k+9]  *  0.0;  )  /*  (1-3,  7-12)  */ 

/*  (xx:  could  do  this  just  at  the  beginning  of  each  trial)  */ 


else  { 


/*  TPBV  guidance;  OPTG  (sans  squib  effects) 


Lcs  -  laiiit4]*c_gamma  -  lamI5]*s_gamma; 

Lsc  ■  lam[4]*s~gaiiuna  -f  lam[Si*c~gamma; 

/*  calculate  commanded  accelerations  */ 
if  (Tau“»0.0)  (  /*  kinematic  eqns  */ 

tmp  «»  V_dot*c_gamma  -  V*gamma_dot*s_gamma; 
tmp2  ■  V*  psi2dot*c_gamma; 

Uctll  ■  tmp*8_psi  +  tmp2*e_p8i; 

Uc[2]  -  tmp*c_p8i  -  tmp2*s_psi; 

DC [3]  ■  V_dot*8  gamma  +  V*gamma_dot*c_gainma; 

) 

else  {  /*  variational  eqns  •/ 

TT  ■  Tau*Tau/2.0; 

FOR3(  det[k]  ■  d[k]  +  Tau*e(kl  +  TT*uTm[k);  )  /*  measured  uT  */ 

/*  (TT,  det  required  only  when  Tau!=0)  */ 

Uc[ll  ■  Lc8*s_psi  +  lam[6)*c_psi;  /*  (partial)  */ 

UcI21  ■  Lc8*c_p8i  -  Iaml61*8_p8i; 

Uc(3]  -  Lsc; 

tmp  -  TT*W[5]*Tau;  /*  0.5*W(5]*Tau*(4-ss)  */ 

FOR3(  Uctkj  +■  W(k]*vlk]  +  dettkl*WI5]*Tau  -  lam[k+9]; 

Uclkj  /■  tmp;  ) 

) 

FOR3(  dbg(  "Uc%d;  %g  ",  k,Uclk)  );  >  dbg("\n"); 

Occs  ■  Uctll*cjpsi  -  Uc[21*s_psi; 

Ocsc  ■  Uctli*8_psi  +  0c[2 l*c_psi; 

L26  >■  m*lam[2]  -f  lam(61: 

1.35  -  m*lamt3]  +  lam[5];  /*  (U4  not  needed)  */ 

lamdotd]  ■  (  rho*V*S*CD0*lamtl]  +  L35*ganuna_dot  )/m; 

lamdot[2]  •=  (  1,26"  (  V"gamma_dot*s_ganuna  -  V_dot*c_gamma  ) 

-  I.C8*Uccs  +  lain[6]*Ocsc  )/(  m*V*c_gamma  ); 

lamdot[3]  ■  (  m*g*(  lamll ] •c_gamma  -  laml 3 ] *8_qamma  ) 

-  L35*V_dot  +  Lsc"Ucsc  -  Lc8*Uc[3]  )/(  m*V  ); 

if  (Tau-«0.0)  FOR3(  lamdot(k+61  «  0.0;  ) 

else  {  tmp  -  2.0*Wt5]/Tau;  /*  (2.0*W[5 ] *Tau*-ss )  */ 

FOR3(  lamdotCk+6]  »  tmp*(  detik]  -  TT*0c[k]  );  ) 

) 

FOR3(  lamdot[k+9]  »  W(kl*Uclk)  +  Tau*lamdot[k+6]  -  lam[k+61  +  SW[k];  ) 

)  /*  (G_mode)  */ 

)  /*  dyn  V 


106 


/*  quad.c 


•include  "standard. h" 

/* 

standard  fns,  consts 

*/ 

/*  Global  variables 

from 

main  */ 

extern  int  debug; 

/* 

debug  flag 

*/ 

extern  int  Abort; 

/* 

fatal  error  flag 

V 

quad(  A,B,C,  proot,mroot,  neg  )  /* 

solve  quadradic 

*/ 

double  A,B,C; 

/* 

input:  Axx+Bx+C»0 

*/ 

double  *proot, *mroot; 

/* 

output:  +,-  roots 

*/ 

int  *neg; 

/* 

output:  error  flag 

*/ 

\ 

double  tmp,  discriin« 

RR; 

/* 

work  vars 

*/ 

/* - 

- */ 

"neg  -  FAI£E; 

/* 

until  proven  otherwise 

*/ 

if  (  A--0.0  ) 

/* 

solution  can  be  linear 

*/ 

if  (  B— 0.0  ) 

Abort 

-  TRUE; 

/*  (stop  trial) 

*/ 

else  *proot  »  ‘mroot  «  -C/B; 

else  (  /*  solution  can  be  quadratic  */ 

RR  *!  -B/(2.0*A);  /*  real  part  (if  discrim  neg )  */ 
discrim  -  RR*RR  -  C/A; 
if  (discrim  <  0.0)  { 

*neg  »  TRUE;  *proot  »  *mroot  *  RR;  /*  use  real  part  only  */ 
dbg(  "\n — >  Warning:  Negative  discriminant:  %g  < — \n",  discrim  ); 
dbg(  "(Components:  a:  %g;  b:  %g;  c:  %g)\n",  A,B,C  ); 
dbg(  "(linear  (-c/b):  %g;  using  real  (-b/2a):  %g)\n",  -C/B,  RR); 

) 

else  {  /*  real  roots  */ 

tmp  ■  8qrt(  discrim  );  *proot  ■  RR  +  tmp; 

"mroot  »  RR  -  tmp; 

)  /*  (discrim)  */ 

)  /*  (A)  V 

dbg(  "plus  root;%g,  minus  root:%g\n",  "proot,  "mroot  ); 

)  /*  quad  */ 
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/•  rng.c  ■/ 

double  gauss (  std  ) 

double  std; 

•define  NUMRNDS  5 


pseudo-gauss ian  pseudo-random  number  gen'r 
Barron  Associates,  Inc.  1/9/87 

desired  gaussian  standard  deviation 
•uniform  random  numbers  to  sum  to  approx 


extern  double  sqrtO; 

double  randonO,  stduse,  sum;  register  int  i; 

stduse  •  std  *  sgrtO.O/MUMRMDS); 

for  (suffi>0.0,  i-0;  i  <  NUMRNDS;  L**) 
sum  random (-stduse,  stduse); 

return (  sum  ); 

)  /*  (gauss)  */ 

double  random(  lo,  hi  )  /*  uniform  pseudo-random  number  generator  */ 

double  lo,  hi;  /*  limits  on  returned  number  (order-invariant)  */ 

{ 

extern  int  srand(),  randO; 
double  tmp,  prop; 
static  int  f irst_time«l ; 

if  (first_time)  {  8rand(12347) ;  f ir8t_time»0;  )  /*  (seed)  */ 

if  (lo  >  hi)  {  tmp«lo;  lo»hi;  hi«tmp;  )  /*  (swap)  */ 

do  prop  ■  ((double)  rand())  /  32767.0;  /*  rand  range:  0-2“15-l  */ 

while  (  (prop  <  0.0)  ||  (prop  >  1.0)  );  /*  xx:  shouldn't  be  nec.*/ 


return(  lo  +  prop*(hi-lo)  ); 


/*  Distribute  between  limits 


1 

iS' 


/*  standard. h  */ 


8 

Si 


§ 


S 

I 

a 


•^5 


tinclude  <stdio.h>  /*  standard  i/o  functions 

/*  (printfO,  fprintfOt  getsO,  fflushO,  strcpyO,  strcatO) 
extern  double  sinO,  cosO,  fabsO,  sqrtO;  /*  math.h  functions 
/*  (also:  tanO,  as  in  ( ) ,  acos  ( ) ,  atan  ( ) ,  exp( ) ,  loglO  ( )  powO,  fnodO) 


•define  PI  3.141S9265359 
•define  DEGtoRAD  (PI/180.) 

•define  rad(h)  DEGtoRAD* (h) 

•define  deg(h)  (h)/DEGtoRAD 

•define  min (a, b)  (  ( (a)< (b) )? (a) : (b)  ) 
•define  max(a,b)  (  ( (a)>(b) )?(a) : (b)  ) 
•define  nl  printf("\n") 

•  define  signum(x)  (((xXO.O)  ?  -1:  1) 
•define  EPS  0.00001 

/*  define  epsig(x)  ( ( f abs (x XEPS )  ?  0 


*/ 

*/ 

V 

*/ 


/* 

degree 

to  radians  converter 

*/ 

/* 

convert 

to  radians 

*/ 

/* 

convert 

to  degrees 

*/ 

/* 

newline 

*/ 

/* 

signum  fn 

with  0  positive 

•/ 

/* 

epsilon  ( 

"small  enough") 

*/ 

s ignum (X  )  ) 

*s ignum  w/  0  band 

•/ 

•define  limit (h, lo, hi )  h  -  ( h< ( lo)? do) : (  h>(hi)?(hi):h  ))  /*  bi-limit*/ 
•define  RSS(a,b)  sqrt ( (a) * (a)+ (b)* (b) >  /*  root  sum  squared  •/ 
•define  RSSS(a,b,c)  sqrt ( (a ) * (a )+ (b)* (b)+ (c ) * (c ) )  /*  RSS  w/  3  */ 


•include  <environ.h> 
•define  SCR_HEIGHT  24 
•define  SCR_WIDTH  80 
•include  <dbg_off.h> 


/*  machine-dependent  params  */ 

/*  debug  macros  (incl.  dt>g )  •/ 
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/*  integrate. h  */ 

typedef  struct  { 

double  u,  up; 
}  STATE_VAR; 

tdefine  N  STATES  22 


/*  a  state  variable 
/*  current  u,  u  dot  vector 


/*  number  of  system  states 


/*  aliases  for  states  and  derivatives  (use  structure  "state") 
tdefine  St(h)  state[h].u 

tdefine  St_DOT(h)  state[h].up 

int  k;  /*  macro  loop  index  */ 

tdefine  L(X)P_St(h)  for  (Jc-0;  )c<N_STATES;  )c++)  (h) 
tdefine  F0R3(h)  for  ()c»l;  k<*3;  k++)  (h) 


V\' 


I 


I 
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/***  system. h  for  129  Tank  projectile  •**/ 

/*  system  constants  •/ 

•define  m  0.966  /*  mass  (slugs)  (31.1  lbs) 

•define  S  0.12151  /*  reference  area  (ft*2) 

/*  (see  also  Mach-related  CDO,  K,  CNa  tables  in  dyn) 

/*  environment:  */ 

•define  g  32.199  /*  gravity  */ 

•define  rho  0.00237692  /*  atmospheric  density  */ 


I  09 

nnuauniM)unu7 


'* 

* 

I 

I 


/*  init.h: 

interface  variables 

*/ 

/*  (common 

variables  in  " 

fandb. h 

")  •/ 

/*  Bacicward 

flight 

*/ 

static  struct 

DESCRIPTION 

backinfoiJ  *  ( 

t include  "fandb.h" 

"N", 

«dM  f. 

200., 

"Final  Isquibs", 

N  n 

9 

)  ; 

NULL, 

0., 

" (end_of_structure  pointer)' 

/*  Forward 

flight 

•/ 

static  struct 

DESCRIPTION 

fwdinfotj  •  { 

< include  "fandb.h" 

"lam4". 

8lam[ 4] , 

200., 

"V_dot  lambda  (const)". 

"phi". 

Sphi, 

0., 

"body  roll  attitude  (deg)”, 

"V", 

SV,  3821.5, 

"velocity  magnitude  (ft/sec)" 

"psi", 

Spsi, 

0., 

"velocity  vector  heading  (deg 

"laml". 

4lam[ 1 ] , 

100., 

"V_dot  lambda". 

•lam2". 

4lam[2] , 

100.  , 

"psi_dot  lambda". 

"lam3", 

4lam(3], 

■100., 

"gamma_dot  lambda". 

"lam7", 

4lam[ 7 ] , 

0., 

"xl_dot  l2unbda". 

"lam8". 

tlam[81 , 

0., 

"x2_dot  lambda". 

"lam9", 

fclamt  9] , 

0., 

"x3_dot  lambda". 

"lamlO", 

6lam[ 10 ] , 

0., 

"vl_dot  lambda". 

"lamll", 

4lam[ll] , 

0., 

"v2_dot  lambda". 

"lainl2", 

4lain[121, 

0., 

■v3_dot  lambda". 

"W6*, 

(W(6], 

0., 

"downrange  weight". 

"W7", 

*.W(7), 

0., 

"crossrange  weight". 

"W8", 

6W(8], 

0., 

"altitude  weight". 

M  M 

9 

NULL, 

0., 

"(end_of_structure  pointer)" 

9 

)  ; 


/*  "use_APN",  Sduse_APN,  0.,  "Use  aPN{s)  (0:no;  l:yes)",  */ 
/*  (XX  make  W6,7,8  bi-directional)  */ 


I 
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fandb. h 


"duration". 

(duration,  10. 

"dt". 

«dtO, 

.05 

•1 

"print_int" 

,  (print  int,  .05 

i 

"integ^. 

(dinteg. 

0., 

"wide" , 

(dwide. 

0., 

"debug". 

(ddebug , 

0., 

"hist" , 

(dhist. 

0., 

1 

"guide", 

(dG_mode, 

2., 

"bounce" , 

(dbounce. 

-1., 

"thresh". 

(Nd  thresh,  0., 

m 

"record". 

(dR_mode, 

0., 

S 

"search". 

(dS~mode, 

0., 

"P", 

(P, 

21., 

& 

"gamma". 

(g5unma. 

2., 

"xTl", 

(xTtll , 

100. , 

"xT2", 

(XT(2] , 

10000., 

"xT3", 

(XT(31, 

200.  , 

"vTl", 

(VT(1]  , 

0., 

"vT2", 

(VT121  , 

0., 

"vT3", 

(VT(31, 

0., 

"uTl", 

(UT(11  , 

0., 

% 

"uT2", 

(UT[2], 

0., 

"uT3", 

(uT[31, 

0., 

"Umean", 

(Umean , 

0., 

"Usig", 

(Usig, 

0., 

"WO  ", 

(W(0] , 

400.  , 

"Wl", 

(W[l], 

0., 

"W2*, 

(W(2), 

0., 

"W3", 

(W(3] , 

0., 

"W4", 

(W[4J, 

100. , 

V 

"W5", 

(W151, 

0.001, 

"lam6". 

(lam[ 6  ] , 

10000. , 

"lamS", 

(lami 5  ] , 

10000., 

"maximum  duration  of  flight  (sec)", 

"timestep  magnitude", 

"printing  interval  (sec)", 

"integration  method  (0:rectangular;  l:Heun)' 
"wide  screen  flag  (0:  80;  1:  132  columns)", 
"debug  flag  (0: normal;  1: debug)", 

"variable  history  flag  (0;none;  l:record)". 


"guidance  (0:ballistic;  Irbaseline;  2:0PTG)", 
"(negrno  bounce;  0-2:bounce  with  that  guidance!' 
"N_dot  threshold  above  which  firing  occurs", 
"recording  (Orend  states;  ltduring;  2:study!", 
"0:none;  l:grid;  2;grope;  3:gauss;  4:unif;  5:fi; 


"roll  rate  (Hz)", 

"velocity  vector  inclination  (deg)". 


"Target  eastward  position  (ft)", 

"Target  northward  position  (ft)", 

"Target  altitude  (ft)", 

"Target  eastward  velocity  (ft/sec)", 

"Target  northward  velocity  (ft/sec)", 

"Target  rate  of  altitude  increase  (ft/sec)", 
"Target  eastward  acceleration  (ft/sec*2)", 
"Target  northward  acceleration  (ft/sec*2)", 
"Target  altitude  acceleration  (ft/sec*2)", 
"mean  target  acceleration  noise  (f/sec*2)", 
"target  acceleration  noise  sigma  (f/sec*2)", 


"N_dot  scaling  weight", 

"Eastward  )cinetic  energy  use  weight", 
"Northward  )cinetic  energy  use  weight", 
"Altitude  Icinetic  energy  use  ./eight", 
"Squib  use  (N)  weight", 

"Predictive  error  penalty  weight". 


"psi_dot  lambda  (const)", 
"gamma_dot  ljunbda  (const)". 


gw 


VVVJ  •cKiwaa 


/•  file:  bounds. h 


/ 


Idef ine 

numF  3 

/*  vary  fwd: 

l<un4 

Iam5 

Iam6 

static 

double 

MinF[numFl  =  { 

-150000. 

,  -150000 

,  -150000 

static 

double 

MaxPinumFl  =  ( 

200000. 

,  100000 

,  150000 

static 

double 

IncF[numFl  =  { 

50000. 

,  50000 

,  50000 

Idefine 

numB  3 

/*  vary  bkwd: 

lam  5 

lam  6 

W4  */ 

static 

double 

MinBInumB]  »  ( 

-10000., 

-10000., 

100.  ); 

static 

double 

MaxBlnumBl  -  { 

10000., 

10000. , 

600.  ); 

static 

double 

IncBinumB]  =  { 

4000., 

4000., 

100.  ); 

I -  129  3<3of  makefile - + 

t  make  -  compile  the  program 

<  make  lint  -  lint  in  the  background 

#  make  go  -  run  the  program 

<  make  all  -  lint  in  the  background  and  then  run 

#  (note:  gives  you  too  many  prompts  for  unknown  reason) 

I - 

fmakefile  template  for  Barron  Associates,  Inc.  7/30/87  (dwa,  ph) 

#n^une  of  final  application 
PROG  >  3dof.out 

tname  of  each  object  file 

OBJ  >  main.o  fly.o  dyn.o  advise. o 

OBJ2  »  rng.o  guad.o  grope. o  step.o  matrix. o 

talso  alt.o 

tlist  of  standard  libraries  to  use 
STDLIBS  =  -Ic  -Im 

Iname  of  BAI  library(ies)  (now  in  standard  place  -  8/3/87/ph) 

IBAILIBS  =  -p  /usr/bai/lib/libBAI.a 
BAILIBS  -  -IBAI 

f location  of  BAI  include  files 
INC  ■  /usr/include 

<  - commands - 

#the  default  (just  type  'make') 

S(PROG):  $(0BJ)  $(0BJ2)  /usr/lib/libBAI.a 

CC  -O  $%  $(0BJ)  $(OBJ2)  $(STDLIBS)  $(BAILIBS) 

#  'make  go'  will  run. 
go:  $(PR0G) 

nice  $ (PROG) 


f  'make  lint'  will  strongly  lint  all  in  background 
lint: 

back  lint  -abchx  *.c 

#  'make  all'  will  lint  then  run 
all:  lint  go 


# - 

-  dependencies  — 

main . o: 

integrate. h 

standard. h 

init. h 

fandb. h 

$(INC)/interfacel.h 

fly.o: 

integrate. h 

standard. h 

system. h 

advise. o: 

standard. h 

bounds . h 

$(INC)/interfacel .h 

dyn .0: 

standard. h 

system. h 

guad.o: 

standard. h 

grope. Or 

step.o: 

standard.h 

grope. h 

matrix. 0: 

grope. h 

f  alt.o: 

f  apnO.o,  apni.o: 


standard. h 


apn.  h 


